数值分析实验报告-插值、三次样条

1、实验报告:顿差值多项式&三次样条1问题：在区间-1,1上分别取n=10、20用两组等距节点对龙格函数f(x) -作多项式插25x值及三次样条插值，对每个n值，分别画出插值函数及f(x)的图形。加深对多项式插值的理解。应用所实验目的：通过编程实现牛顿插值方法和三次样条方法, 编程序解决实际算例。实验要求：认真分析问题，深刻理解相关理论知识并能熟练应用; 编写相关程序并进行实验； 调试程序，得到最终结果；分析解释实验结果； 按照要求完成实验报告。实验原理：详见数值分析 第5版第二章相关内容。实验内容：(1 )牛顿插值多项式1.1 当 n=10 时：在Matlab下编写代码完成计算和画图。结

2、果如下: 代码：clear allclc x1=-1:0.2:1;y1=1./(1+25.*x1A2);n=len gth(x1);f=y1(：);for j=2: nfor i=n :-1:jf(i)=(f(i)-f(i-1)/(x1(i)-x1(i-j+1);end endsyms F X p;F(1)=1; p(1)=y1(1);for i=2: nF(i)=F(i-1)*(x-x1(i-1);P (i)=f(i)*F(i);endsyms PP=sum (p);P10=v pa(ex pan d( P),5); xO=-1:O.OO1:1;y0=subs( P,x,xO); y2=su

3、bs(1/(1+25*xA2),x,xO);p lot(xO,yO,xO,y2) grid on xlabel( ylabel('x')'y')即我们所求的牛顿插值多项式，其结果为：P1O( x)=-22O.94*x1O+494.91*x8-9.5O65e-14*x7-381.43*x6-8.5O4e-14*x5+123.36*x4+2.O2O2e-1P104*xA3-16.855*xA2-6.6594e-16*x+1.0并且这里也能得到该牛顿插值多项式的在-1,1上的图形，并和原函数进行对比（见Fig.1)。A111 IIID£口吕Fig.1牛顿插值

4、多项式（n=1O）函数和原函数图形 从图形中我们可以明显的观察出插值函数在两端点处发生了剧烈的波动， 误差，即龙格现象，当n=2O时，这一现象将更加明显。1.2 当 n=2O 时： 对n=1O的代码进行修改就可以得到n=2O时的代码。将"x1=-1:O.2:1; ”改为“x1=-1:O.1:1;”即可。运行程序，我们得到n=2O时的牛顿插值多项式，结果为：P2O（x）=26O188.O*x2O - 1.O121e6*x18 + 2.6193e-12*x17 + 1.6392e6*x16 + 2.248e-11*x15产生了极大的-1.4429e6*x14 - 4.6331e-11*x

5、13+ 757299.0*x12 + 1.7687e-11*x11- 245255.0*x10+ 2.1019e-11*xA9 + 49318.0*xA8 + 3.5903e-12*x7 - 6119.2*x6 - 1.5935e-12*x5 +470.85*x4 + 1.3597e-14*x3 - 24.143*x2 - 1.738e-14*x + 1.0同样的，这里得到了该牛顿插值多项式的在-1，1上的图形，并和原函数进行对比（见Fig.2)。Fig.2牛顿插值多项式（n=20）函数和原函数图形当n=20时，端点处发生了更加剧烈的震荡。表明随着分段不断增加，牛顿插值多项式 与原函数的误差不

6、但没有减少，反而变得更大了。（2 ）三次样条2.1 当 n =10 时：在Matlab下编写代码完成计算和画图。结果如下: 代码：clear all cic x1=-1:0.2:1;y1=1./(1+25.*x1A2);syms x m1=subs(diff(1/(1+25*x2),-1);m2=subs(diff(1/(1+25*xA2),1);n=len gth(x1);syms a b h f d for i=1: n-1h(i)=x1(i+1)-x1(i);f(i)=(y1(i+1)-y1(i)/(x1(i+1)-x1(i);end a(n)=1; b(1)=1;for i=2:n-1

7、a(i)=h(i-1)/(h(i-1)+h(i); b(i)=h(i)/(h(i-1)+h(i);end d(1)=6/h*(f(1)-m1); d(n )=6/h( n-1)*(m2-f( n-1);for i=2:n-1d(i)=6*(f(i)-f(i-1)/(h(i-1)+h(i);endD=d：A=2.*eye( n); for i=1: n-1A(i,i+1)=b(i);A(i+1,i)=a(i+1);endM=AA-1*D; for i=1: n-1 s(i)=M(i)*(x1(i+1)-x)A3/h(i)/6+M(i+1)*(x-x1(i)A3/h(i)/6+(y1(i)-M(i

8、)* h(i)A2/6)*(x1(i+1)-x)/h(i)+(y1(i+1)-M(i+1)*h(i)A2/6)*(x-x1(i)/h(i);endS=vp a(ex pan d(s.'),5);for i=1: n-1x0=-1- (2/( n-1)+(2/( n-1)*i:0.001:-1+(2/( n-1)*i; y0=subs(s(i),x,x0);Fig.5)。plot(x0,y0) holdon-1,-0.8-0.8,-0.6-0.6,-0.4-0.4,-0.2-0.2,00,0.20.2,0.40.4,0.60.6,0.80.8,1-1, 1上的图形，并和原函数进行对比(见

9、end y2=subs(1/(1+25*x2),x,-1:0.001:1);plot(-1:0.001:1,y2,'r')grid on xlabel('x')ylabel('y')S即为我们所求的三次样条，其结果为：So(x)=0.08225*x3+0.36953*x2+0.56627*x+0.317450.96279*xA3+2.4828*xA2+2.2569*x+0.76829 0.81773*xA3+2.2217*xA2+2.1002*x+0.7369613.413*x3+17.336*x2+8.1461*x+1.5431-54.471*

10、x3-23.394*x2-1.8741e-17*x+1.054.471*x3-23.394*x2+1.9683e-17*x+1.0 -13.413*x3+17.336*x2-8.1461*x+1.5431-0.81773*x3+2.2217*x2-2.1002*x+0.73696 -0.96279*x3+2.4828*x2-2.2569*x+0.76829-0.08225*x3+0.36953*x2-0.56627*x+0.31745并且这里可以得到该三次样条的在Fig.3三次样条（n=10）函数和原函数图形 从图形我们可以看出，三次样条图形和原函数图形非常接近，误差相对较小。2.2 当 n=

11、20 时：同样的，将上面代码中的“ x1=-1:0.2:1; 我们可以得到n=20时的三次样条，其结果为”改为“ x1=-1:0.1:1;”即可。运行程序,S2o(x)=0.16461*xA3+0.59746*xA2+0.77505*x+0.380660.27191*xA3+0.88717*xA2+1.0358*x+0.458880.46379*x3+1.3477*x2+1.4042*x+0.55712-1,-0.9-0.9,-0.8-0.8,-0.70.86962*x3+2.1999*x2+2.0008*x+0.696321.5804*x3+3.4792*x2+2.7683*x+0.8498

12、43.5442*x3+6.425*x2+4.2412*x+1.09535.7284*x3+9.046*x2+5.2896*x+1.235112.534*x3+15.171*x2+7.1273*x+1.4189-0.7,-0.6-0.6,-0.5-0.5,-0.4-0.4,-0.3-03-0.2-32.789*x3-12.023*x2+1.6884*x+1.0563-89.07*x3-28.907*x2+5.1067e-17*x+1.0-02-0.1-0.1,089.07*x3-28.907*x2-3.9148e-17*x+1.032.789*x3-12.023*x2-1.6884*x+1.0563-12.534*x3+15.171*x2-7.1273*x+1.41890,0.10.1,0.20.2,0.3-5.7284*x3+9.046*x2-5.2896*x+1.2351-3.5442*x3+6.425*x2-4.2412*x+1.09530.3,0.40.4,0.5-1.5804*x3+3.4792*x2-2.7683*x+0.84984-0.86962*x3+2.1999*x2-2.0008*x+0.696320.5,0.60.6,0.7-0.46379*x3+1.3477*x2-1.4042*x+0.55712-0.271