勒让德逼近编程实现“勒让德逼近”.doc

CENTRAL SOUTH UNIVERSITY MATLAB程序设计实践报告 勒让德逼近编程实现“勒让德逼近”，并用勒让德公式（取6项）逼近函数：1/(2-x)，并求出当x=0.5时的函数值。 算法说明： 当一个连续函数定义在区间[-1，1]上时，它可以展开成勒让德级数。即： 其中Pn(x)为n次勒让德多项式 勒让德的逼近要求被逼近函数定义在区间【-1,1】上，勒让德多项式也可以通过递推来定义： 它们之间满足如下的正交关系： 在实际应用中，可以根据所需的精度来截取有限项数。勒让德级数中的系数由下式确定： 在MATLAB中编程实现的勒让德逼近法函数为：Legendre。 功能：用勒让德多项式逼近已知函数。 调用格式：f= Legendre（y，k）或f= Legendre（y，k，x0） 其中，y为已知函数； k为逼近已知函数所需项数； x0为逼近点的x坐标： f为求得的勒让德逼近多项式或在x0处的逼近值。 源程序： 第一个源程序： function f=Legendre(y,k,x0) %用勒让德多项式逼近已知函数 %已知函数：y %逼近已知函数所需项数：k %逼近点的x坐标：x0 %求得的勒让德逼近多项式或在x0处的逼近值：f syms t; P(1:k+1)=t; P(1)=1; P(2)=t; c(1:k+1)=0.0; c(1)=int(subs(y,findsym(sym(y)),sym(t))*P(1),t,-1,1)/2; c(2)=int(subs(y,findsym(sym(y)),sym(t))*P(2),t,-1,1)/2; f=c(1)+c(2)*t; for i=3:k+1 P(i)=((2*i-3)*P(i-1)*t-(i-2)*P(i-2))/(i-1); c(i)=int(subs(y,findsym(sym(y)),t)*P(i),t,-1,1)/2; f=f+c(i)*P(i); if(i==k+1) if(nargin==3) %输入参数个数为三 f=subs(f,t,x0); else f=vpa(f,6); %对逼近值f取六位有效数字 end end end 第二个源程序： f=Legendre(1/(2-x),6) f=Legendre(1/(2-x),6，0.5) ③运行结果： f = 0.0955158*t - 0.0000116697*t*(2.2*t*(2.25*t*(4.66667*t - 2.33333*t*(7.5*t^2 - 2.5)) + 10.125*t^2 - 3.375) - 5.86667*t + 2.93333*t*(7.5*t^2 - 2.5)) - 0.0000565953*t*(2.25*t*(4.66667*t - 2.33333*t*(7.5*t^2 - 2.5)) + 10.125*t^2 - 3.375) + 0t*(7.5*t^2 - 2.5) - 0.000275682*t*(4.66667*t - 2.33333*t*(7.5*t^2 - 2.5)) + 0.0305351*t^2 + 0.539128 f = 0.5935 从逼近结果上来看，函数的准确值为1/(2-0.5)=0.6667. ④流程图如下： 第一题流程图： 是 否 是 否 是 第二题流程图： 二．编程解决以下科学计算和工程实际问题。 简支梁受左半均匀分布载荷q及右边L/4处集中力偶M0作用（如下图1-1），求其弯矩、转角和挠度。设L=2m，q=1000N/m，M0=900N*m，E=200*109N/m2,I=2*10-6m4. 图1-1 解题思路： 首先对简支梁进行受力分析，受力分析图（如下图1-2）所示： 图1-2 从材料力学的知识可知道，由弯矩求转角要经过一次不定积分，而由转角求挠度又要经过一次不定积分，通常这是很麻烦而且容易出错的，而在MATLAB中，可用cumsum函数或cumtrapz函数作近似的不定积分，只要x取得足够密，其结果将相当准确，而程序非常简单。 第一步：计算支反力 设支座a和b处的支反力分别为Na和Nb，则据∑Ma=0，∑Fy=0得到平衡方程为： Nb=(q*L^2/8+M0)/L Na=q*L/2-Nb 第二步：建立弯矩方程 以截面c，d为分界面，将梁划分为ac，cd，db三段 分别建立ac，cd，db三段对应的弯矩方程： M1=Na*x-q*x.^2/2; 0≦x≦L/2 M2=Nb*(L-x)