数值分析课程设计

1、【精品文档】如有侵权，请联系网站删除，仅供学习与交流数值分析课程设计.精品文档.2013/2014第一学期数值分析课程设计设计题目： 插值算法总结与Matlab编程实现 信息与计算科学 专业 学号 112086202 姓名 胡波 学号 112086237 姓名 许世开 学号 112086235 姓名 郭海波 成绩 指导老师 郭尊光 摘 要本论文首先分析用插值多项式表示复杂函数f（x）的必要性，然后运用插值的方法根据给定的函数表构造出一个既能反映函数f (x)的特征，又便于计算的简单函数p(x)，使得p(x)近似f (x)。文章通过对拉格朗日插值法存在的缺点的分析，引入了牛顿插值法、埃尔米特插值

2、法、分段低次插值法及三次样条插值法，给出了每种插值算法的算法结构、具体的程序MATLAB代码及误差估计的方法，并且对每种插值方法的优缺点进行了分析。关键词：拉格朗日插值，牛顿插值，埃尔米特插值，分段低次插值，三次样条插值问题提出 设f (x)的函数表为x =2,4,6,8,10;y =9.8,9.2,8.1,6.4,3.8，y=0.96,0.87,0.58,0.45,0.25,g（x）=，(-5x5)。请根据给出的函数及函数表解下面的问题。 (1)证明插值多项式的存在唯一性，以x ,y为节点计算其插值多项式。（2）用基函数构造拉格朗日插值多项式，以x ,y为节点做拉格朗日插值曲线。（3）分析牛

3、顿插值多项式的算法及余项，利用牛顿插值多项式，求f (6.596)并估计误差。 (4）埃尔米特插值多项式的构造及余项，利用埃尔米特插值多项式，求f (6.596)。（5）在g（x）=，(-5x5)的条件下，分析分段线性插值的必要性（龙格现象）及误差。（6）在g（x）=，(-5x5)的条件下，建立三次样条函数，通过观察图象，得出结论“三次样条插值s(x)较十次拉格朗日插值L(x)能更好地逼近被插函数f (x)”。1 插值多项式1.1插值多项式定义设函数y=f (x)在区间a,b上有定义，且已知在点上的值，若存在一个n次的多项式p(x)，使得 （1.1） 成立，则称p(x)为f (x)的n次拉格朗

4、日插值多项式，点称为插值节点，包含插值节点的区间a,b称为插值区间。注： 1、上述定义中的(1.1)式通常称为插值条件，它是n+1个条件。 2、插值多项式的几何意义如下 求曲线y=p(x)，使其通过给定的n+1个点() ( i=0,1, 2,n )，并用它近似已知曲线y=f (x)。1.2插值多项式存在唯一性与计算首先我们证明插值多项式存在且唯一。由插值条件(1.1),有 + + = （1.2）式（1.2）是一个关于未知量的线性方程组，其系数矩阵的行列式是一个范德蒙行列式由于插值节点是互异的，所以V0。从而线性方程组(1.2)有唯一解，p(x)存在且唯一。以上的证明表明：由n+1个插值条件可唯

5、一地确定一个n个拉氏插值多项式。下面我们在问题1的条件下计算以x=2,4,6,8,10,y=9.8,9.2,8.1,6.4,3.8为插值节点的插值多项式。MATLAB代码如下：format shortx=2,4,6,8,10'y=9.8,9.2,8.1,6.4,3.8'A=x.0,x.1,x.2,x.3B=Ay运行结果为：B =10.1600,-0.1476,-0.0071,-0.0042' 这表明: 。2 用基函数构造拉格朗日插值多项式2.1 拉格朗日插值多项式的构造利用求解线性方程组(1.2)求拉格朗日插值多项式不便计算机实现，不易编程。在实际应用中，一般采用基函数

6、法来构造。下面，我们介绍这一方法，这里，我们将n次拉格朗日插值多项式用来记。设,其中是待定的n次多项式。首先构造个插值节点上的插值基函数，对任一点所对应的插值基函数，由于在所有取零值，因此有因子。又因是一个次数不超过的多项式，所以只可能相差一个常数因子，固可表示成：利用得：于是因此满足 的插值多项式可表示为：从而次拉格朗日插值多项式为：2.2拉格朗日插值多项式在MATLAB中的实现及使用在MATLAB中构造拉格朗日插值函数代码如下：function yy=lagelangrichazhi(x,y,xx)m=length(x);n=length(y);if m=n;end s=0; for i=

7、1:n t=ones(1,length(xx); for j=1:n if j=i,t=t.*(xx-x(j)/(x(i)-x(j); end end s=s+t*y(i); end yy=s;利用该函数做题一中x，y的插值曲线代码如下：clearx=2,4,6,8,10;y=9.8,9.2,8.1,6.4,3.8;xi=2:0.1:10;yi=lagelangrichazhi(x,y,xi);plot(x,y,'o',xi,yi,'k');title('拉格朗日插值')得到拉格朗日插值曲线如图2.1所示： 图2.1 拉格朗日插值曲线3牛顿多项式

8、lagrange插值多项式作为一种计算方案，公式简洁，做理论分析也方便。其缺点是，当节点改变时，公式需要重建，计算量大；如果还要根据精度要求，选取合适的节点和插值多项式的次数，则只好逐次计算出 , 等，并做误差试算，才可以做到，这当然是不理想的。为此，人们从改进插值多项式的形式入手，给出另一种风格的插值公式，这就是Newton(牛顿)插值公式。3.1均差概念及其性质设函数在互异节点上的值为 , 等，定义：(1) 在上的1阶均差为 (2) 在上的2阶均差为 （3）递推地，在上的阶均差为同时规定在上的零阶均差为 。均差具有的性质如下：性质1 阶均差可以表示成个函数值的线性组合，即或记为 性质2 均

9、差对其节点是对称的，即节点按任意顺序排列，其值不变。如这个性质称为对称性。性质3 如果是的次多项式，则其1阶均差是的次多项式，且由此递推可得阶均差为零阶均差，阶均差为零。差商表的构造如表1.3所示：表3.1 差商表 零阶均差1阶均差2阶均差 3阶均差 利用MATLAB计算题1中x，y的各阶均差代码如下：x=2,4,6,8,10'y=9.8,9.2,8.1,6.4,3.8'n=max(size(x); %计算节点个数z=zeros(n,n+1);%建一个n 行n+1 列的零矩阵z(:,1)=x; %将x 赋给z 的第1 列z(:,2)=y;%将y 赋给z 的第2 列for j=1

10、:n-1for i=j+1:nz(i,j+2)=(y(i)-y(i-1)/(x(i)-x(i-j); %计算j 阶差商endfor i=j:ny(i)=z(i,j+2); %将j 阶差商值赋给向量y 的对应分量endendz计算结果为：3.2牛顿插值多项式的算法及余项由均差的定义可知对上述第二式两端乘以，第三式两端乘以, 最后一式两端乘以，然后由后一式的左端代入前一式的右端(第二项)，即可得 显然，是次多项式，又有，故有，从而满足条件，这就说明是关于的次插值多项式，通常称它为Newton插值公式，为插值余项。在MATLAB中使用牛顿插值法求得题1中以x，y为插值节点的条件下f(6.596)的值

11、并估计误差的代码如下：format longx=2,4,6,8,10;y=9.8,9.2,8.1,6.4,3.8;n=max(size(x);%计算差商for j=1:1:n-1for i=n:-1:j+1y(i)=(y(i)-y(i-1)/(x(i)-x(i-j);endend%计算牛顿插值多项式的值z=6.596;%需求其函数值的点p=y(n);for i=n-1:-1:1p=y(i)+p*(z-x(i);endfprintf('f(%8.5f)=%8.5fn',z,p);%误差估计r=y(n);for i=1:nr=r*(z-x(i);endr=abs(r);fprint

12、f('截断误差R(%8.5f)=%15.14fn',z,r);计算结果如下：4 埃尔米特插值多项式 拉氏插值多项式的还存在着一个缺陷：在插值节点处，被插函数曲线与插值多项式曲线的切线方向可能会不一致。不少实际的插值问题不但要求在节点上函数值相等，而且还要求对应的导数值也相等，甚至要求高阶导数也相等，满足这种要求的插值多项式称为埃尔米特（Hermite）插值多项式。4.1埃尔米特插值多项式定义设在节点上，若存在一个2n+1次多项式H(x)，使得 成立，则称H(x)为f (x)的2n+1次埃尔米特插值多项式。4.2埃尔米特插值多项式的构造及余项设计2n+2个插值基函数，每一个基函数

13、都是2n+1次多项式，并且满足以下条件于是满足条件(2.5)的插值多项式H(x)可以表示成为下面的问题就是求满足条件(2.6)的基函数。a a 利用拉氏插值基函数以及条件(2.6)，可知的形式应为又整理得解出由于于是类似地，可推导出。仿照拉格朗日余项的证明方法，若f (x)在(a,b)内的2n+2阶导数存在，则其插值余项为在MATLAB中构造埃尔米特函数代码如下： function f=Hermite(x,y,y_1,x0) syms t; f=0.0; if(length(x)=length(y) if(length(y)=length(y_1) n=length(x); else disp

14、('y和y的导数的维数不相等！'); return; end else disp('y和y的导数的维数不相等！'); return; end for i=1:n h=1.0; a=0.0; for j=1:n if(j=i) h=h*(t-x(j)2/(x(i)-x(j)2; a=a+1/(x(i)-x(j); end end f=f+h*(x(i)-t)*(2*a*y(i)-y_1(i)+y(i); if(i=n) if(nargin=4) f=subs(f,'t',x0); else f=vpa(f,6); end end end在MATLA

15、B中计算以x，y为插值节点，且对应的一阶导数为y的埃尔米特多项式，并计算在x=6.596点处函数的值代码如下：clearx=2,4,6,8,10;y=9.8,9.2,8.1,6.4,3.8;y1=0.96,0.87,0.58,0.45,0.25;f=Hermite(x,y,y1)f1=Hermite(x,y,y1,6.596)计算结果为：f=0.000108507*(8.53667*t - 24.9467)*(t - 6.0)2*(t - 10.0)2*(t - 2.0)2*(t - 8.0)2 - 0.000108507*(4.88333*t - 45.4667)*(t - 6.0)2*(t

16、 - 10.0)2*(t - 2.0)2*(t - 4.0)2 + 0.00000678168*(21.3767*t - 32.9533)*(t - 6.0)2*(t - 10.0)2*(t - 4.0)2*(t - 8.0)2 - 0.00000678168*(7.66667*t - 80.4667)*(t - 6.0)2*(t - 2.0)2*(t - 4.0)2*(t - 8.0)2 + 0.000244141*(t - 10.0)2*(t - 2.0)2*(t - 4.0)2*(t - 8.0)2*(0.58*t + 4.62)f1 = 8.19065 分段线性插值和分段三次hermi

17、te插值5.1龙格现象 lagrange插值方法使根据区间a,b上给出的节点构造插值多项式近似f(x)。一般认为的次数n越高，逼近f(x)的效果越好，龙格现象表明 时，不一定收敛到f(x)。我们以函数为例研究龙格现象，MATLAB代码如下：x=-5:1:5.1;y=1./(1+x.2);plot(x,y,'r');n=max(size(x);hold onz=-5:0.005:5; m=max(size(z); for k=1:mL(k)=0;for i=1:nt=1;for j=1:n if j=it=t*(z(k)-x(j)/(x(i)-x(j);endendL(k)=L(

18、k)+t*y(i);endendplot(z,L,'b')hold offtitle('龙格现象 ')图5.1 龙格现象这说明高次插值多项式近似的效果并不好，因而通常并不采用高次插值，而采用分段低次插值。5.2分段线性插值分段线性插值就是通过插值点用折线段连接起来逼近。设已知节点a=x0 <x1 <<xn =b上的函数值yk =f(xk)(k=0,1,n)，已知函数y= f(x) 在这些插值结点的函数值为yk =f(xk)(k=0,1,n)求一个分段函数Ih(x), 使其满足：  (1) Ih(xk )=yk ，(k=0,1

19、,n);  (2) 在每个区间xk ,xk+1 上,Ih (x)是个一次函数。易知,Ih(x)是个折线函数, 在每个区间xk ,xk+1 上,(k=0,1,n) 于是, Ih (x)在a,b上是连续的，但其一阶导数是不连续的。分段线性函数的基函数就是一阶lagrange插值多项式。其误差为（ 其中）在matalb中一为插值函数interp1默认的就是线性插值。在MATLAB中做函数g（x）=，(-5x5)插值节点数不同的插值曲线，可以看到分段低次插值不会产生龙格现象，代码如下：clear%输入插值区间的等分数disp('给出插值区间的等分数n(5 的倍数)'

20、)n=input('n = ');%作被插函数图象u=-5:(10/200):5;v=2*u./(1+u.2);plot(u,v,'r');hold on%固化图形屏幕%给出插值条件x=-5:(10/n):5;y=2*x./(1+x.2);%取需用分段线性插值函数计算其函数值的点zz=-5:10/(2*n):5;m=max(size(z);%计算分段线性插值函数在z处的值，并作图for k=1:mi=1;while i<=nif z(k)>=x(i) && z(k)<=x(i+1)Ih(k)=y(i)*(z(k)-x(i+1)/

21、(x(i)-x(i+1)+y(i+1)*(z(k)-x(i)/(x(i+1)-x(i);breakelsei=i+1;endendendplot(z,Ih,'b')hold off%释放固化的图形屏幕输入等分数为5时，拉格朗日插值曲线如图5.2所示，图5.2 等分数为5时的拉格朗日插值曲线输入等分数为800时，拉格朗日插值曲线如图5.3所示， 图5.3 等分数为800时的拉格朗日插值曲线6 三次样条插值多项式通过平面上在n+1个不同的已知点()，可以作一条光滑的n次代数曲线，但当点很多时，除了插值多项式的次数高，计算复杂外，从逼近的效果来看，随着插值多项式的次数增高，逼近的效果

22、并不是随之越来越理想，相反还会产生龙格现象。采用分段线性插值函数去逼近被插函数，能保持较好的逼近效果，但其曲线在节点处仅仅只能保证连续，却并不是光滑的。在生产和科学实验中，对所构造的插值曲线，既要求简单，又要求在连接处比较光滑。即所构造的分段插值函数在分段上要求多项式次数低，而在分点（节点）上还具有连续的低阶导数，满足这样条件的插值函数，我们称为样条插值。6.1三次样条插值多项式算法组织所谓三次样条插值多项式是一种分段函数，它在节点分成的每个小区间上是3次多项式，其在此区间上的表达式如下：因此，只要确定了的值，就确定了整个表达式，的计算方法如下：令：则满足如下n-1个方程：对于第一种边界条件下

23、有如果令那么解就可以为 6. 2三次样条插值多项式算例分析做以x，y为插值节点的三次样条插值曲线代码如下：x = 2,4,6,8,10;y = 9.8,9.2,8.1,6.4,3.8;xx=2:0.2:10;yy = spline(x,y,xx);plot(x,y,'o',xx,yy)title('三次样条插值曲线')grid on所得结果如图6.1所示：图6.1 三次样条插值曲线6.3三次样条插值函数与拉格朗日插值多项式的效果比较为显示三次样条插值函数与拉格朗日插值多项式的效果比较，我们用MATLAB分别做g（x）=，(-5x5)的函数曲线、拉格朗日插值曲线、

24、三次样条插值曲线，代码如下：clearformat short%作y=2*x/(1+x2)的图象u=-5:(10/200):5;v=2*u./(1+u.2);plot(u,v,'r')%红线表示函数图象hold on%给出插值条件x=-5:1:5;y=2*x./(1+x.2);y_bar(1)=152/(26)2; %在-5,5的一阶导数值y_bar(2)=152/(26)2;%计算每个小区间长度for i=1:10h(i)=x(i+1)-x(i);endfor i=2:10lmd(i)=h(i)/(h(i-1)+h(i);mu(i)=1-lmd(i);b(i)=6*(y(i+

25、1)-y(i)/h(i)-(y(i)-y(i-1)/h(i-1)/(h(i)+h(i-1);endlmd(1)=1;mu(11)=1;b(1)=(6/h(1)*(y(2)-y(1)/h(1)-y_bar(1);b(11)=(6/h(10)*(y_bar(2)-(y(11)-y(10)/h(10);%用追赶法解三对角方程,求矩mbata(1)=lmd(1)/2;for i=2:10bata(i)=lmd(i)/(2-mu(i)*bata(i-1);endu(1)=b(1)/2;for i=2:11u(i)=(b(i)-mu(i)*u(i-1)/(2-mu(i)*bata(i-1);endm(11

26、)=u(11);for i=10:-1:1m(i)=u(i)-bata(i)*m(i+1);end%给出需计算函数值的点z=-5:0.5:5;n=max(size(z);%计算样条函数的值for j=1:ni=1;while i<=10if z(j)>=x(i) & z(j)<=x(i+1)s(j)=m(i)*(x(i+1)-z(j)3/(6*h(i)+m(i+1)*(z(j)-x(i)3/(6*h(i)+(y(i)-m(i)*h(i).2/6)*(x(i+1)-z(j)/h(i)+(y(i+1)-m(i+1)*h(i)2/6)*(z(j)-x(i)/h(i);bre

27、akelsei=i+1;endendend%计算拉氏插值多项式在z处的函数值for k=1:nL(k)=0;for i=1:11t(i)=1;for j=1:11if j=it(i)=t(i)*(z(k)-x(j)/(x(i)-x(j);endendL(k)=L(k)+t(i)*y(i);endend%计算被插函数在z处的值q=2*z./(1+z.2);%输入值与图象,作比较fprintf(' 节点 函数值 样条值 拉氏插值n')disp(z',q',s',L')plot(z,L,'b')%蓝线表示拉氏插值图象plot(z,s,&

28、#39;y')%黄线表示三次插值样条图象hold off所得结果如图6.2所示图6.2 三次插值样条与十次拉格朗日插值效果比较输出结果如下： 节点 函数值 样条值 拉氏插值 -5.0000 -0.3846 -0.3846 -0.3846 -4.5000 -0.4235 -0.3757 0.2572 -4.0000 -0.4706 -0.4706 -0.4706 -3.5000 -0.5283 -0.5455 -0.7007 -3.0000 -0.6000 -0.6000 -0.6000 -2.5000 -0.6897 -0.6682 -0.5970 -2.0000 -0.8000 -0.8000 -0.8000 -1.5000 -0.9231 -0.9905 -1.0195 -1.0000 -1.0000 -1.0000 -1.0000 -0.5000 -0.8000 -0.6198 -0.6264 0 0 0 0 0.5000