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

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

2、llclcx1=-1:0.2:1;y1=1./(1+25.*x1.2);n=length(x1);f=y1(:);for j=2:n for i=n:-1:j f(i)=(f(i)-f(i-1)/(x1(i)-x1(i-j+1); endendsyms F x p;F(1)=1;p(1)=y1(1);for i=2:n F(i)=F(i-1)*(x-x1(i-1); p(i)=f(i)*F(i);endsyms PP=sum(p);P10=vpa(expand(P),5);x0=-1:0.001:1;y0=subs(P,x,x0);y2=subs(1/(1+25*x2),x,x0);plot(

3、x0,y0,x0,y2)grid onxlabel(x)ylabel(y)P10即我们所求的牛顿插值多项式，其结果为：P10（x）=-220.94*x10+494.91*x8-9.5065e-14*x7-381.43*x6-8.504e-14*x5+123.36*x4+2.0202e-14*x3-16.855*x2-6.6594e-16*x+1.0并且这里也能得到该牛顿插值多项式的在-1，1上的图形，并和原函数进行对比（见Fig.1）。Fig.1 牛顿插值多项式（n=10）函数和原函数图形从图形中我们可以明显的观察出插值函数在两端点处发生了剧烈的波动，产生了极大的误差，即龙格现象，当n=20时

4、，这一现象将更加明显。1.2 当n=20时：对n=10的代码进行修改就可以得到n=20时的代码。将“x1=-1:0.2:1;”改为“x1=-1:0.1:1;”即可。运行程序，我们得到n=20时的牛顿插值多项式，结果为：P20(x)= 260188.0*x20 - 1.0121e6*x18 + 2.6193e-12*x17 + 1.6392e6*x16 + 2.248e-11*x15 - 1.4429e6*x14 - 4.6331e-11*x13 + 757299.0*x12 + 1.7687e-11*x11 - 245255.0*x10 + 2.1019e-11*x9 + 49318.0*x8

5、 + 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时，端点处发生了更加剧烈的震荡。表明随着分段不断增加，牛顿插值多项式与原函数的误差不但没有减少，反而变得更大了。（2）三次样条2.1 当n=10时：在Matlab下编写代码完成计算和画图。结果如下：代码：clear allclcx1=-1

6、:0.2:1;y1=1./(1+25.*x1.2);syms xm1=subs(diff(1/(1+25*x2),-1);m2=subs(diff(1/(1+25*x2),1);n=length(x1);syms a b h f dfor i=1:n-1 h(i)=x1(i+1)-x1(i); f(i)=(y1(i+1)-y1(i)/(x1(i+1)-x1(i);enda(n)=1;b(1)=1;for i=2:n-1 a(i)=h(i-1)/(h(i-1)+h(i); b(i)=h(i)/(h(i-1)+h(i);endd(1)=6/h(1)*(f(1)-m1);d(n)=6/h(n-1)*

7、(m2-f(n-1);for i=2:n-1 d(i)=6*(f(i)-f(i-1)/(h(i-1)+h(i);endD=d;A=2.*eye(n);for i=1:n-1 A(i,i+1)=b(i); A(i+1,i)=a(i+1);endM=A-1*D;for i=1:n-1 s(i)=M(i)*(x1(i+1)-x)3/h(i)/6+M(i+1)*(x-x1(i)3/h(i)/6+(y1(i)-M(i)*h(i)2/6)*(x1(i+1)-x)/h(i)+(y1(i+1)-M(i+1)*h(i)2/6)*(x-x1(i)/h(i);endS=vpa(expand(s.),5);for i

8、=1:n-1 x0=-1-(2/(n-1)+(2/(n-1)*i:0.001:-1+(2/(n-1)*i; y0=subs(s(i),x,x0); plot(x0,y0) hold onendy2=subs(1/(1+25*x2),x,-1:0.001:1);plot(-1:0.001:1,y2,r)grid onxlabel(x)ylabel(y)S即为我们所求的三次样条，其结果为：S10(x) =0.08225*x3+0.36953*x2+0.56627*x+0.31745-1,-0.80.96279*x3+2.4828*x2+2.2569*x+0.76829-0.8,-0.60.8177

9、3*x3+2.2217*x2+2.1002*x+0.73696-0.6,-0.413.413*x3+17.336*x2+8.1461*x+1.5431-0.4,-0.2-54.471*x3-23.394*x2-1.8741e-17*x+1.0-0.2,054.471*x3-23.394*x2+1.9683e-17*x+1.00,0.2-13.413*x3+17.336*x2-8.1461*x+1.54310.2,0.4-0.81773*x3+2.2217*x2-2.1002*x+0.736960.4,0.6-0.96279*x3+2.4828*x2-2.2569*x+0.768290.6,0.

10、8-0.08225*x3+0.36953*x2-0.56627*x+0.317450.8,1并且这里可以得到该三次样条的在-1，1上的图形，并和原函数进行对比（见Fig.3）。Fig.3 三次样条（n=10）函数和原函数图形从图形我们可以看出，三次样条图形和原函数图形非常接近，误差相对较小。2.2 当n=20时：同样的，将上面代码中的“x1=-1:0.2:1;”改为“x1=-1:0.1:1;”即可。运行程序，我们可以得到n=20时的三次样条，其结果为S20(x) =0.16461*x3+0.59746*x2+0.77505*x+0.38066-1,-0.90.27191*x3+0.88717*

11、x2+1.0358*x+0.45888-0.9,-0.80.46379*x3+1.3477*x2+1.4042*x+0.55712-0.8,-0.70.86962*x3+2.1999*x2+2.0008*x+0.69632-0.7,-0.61.5804*x3+3.4792*x2+2.7683*x+0.84984-0.6,-0.53.5442*x3+6.425*x2+4.2412*x+1.0953-0.5,-0.45.7284*x3+9.046*x2+5.2896*x+1.2351-0.4,-0.312.534*x3+15.171*x2+7.1273*x+1.4189-0.3,-0.2-32.7

12、89*x3-12.023*x2+1.6884*x+1.0563-0.2,-0.1-89.07*x3-28.907*x2+5.1067e-17*x+1.0-0.1,089.07*x3-28.907*x2-3.9148e-17*x+1.00,0.132.789*x3-12.023*x2-1.6884*x+1.05630.1,0.2-12.534*x3+15.171*x2-7.1273*x+1.41890.2,0.3-5.7284*x3+9.046*x2-5.2896*x+1.23510.3,0.4-3.5442*x3+6.425*x2-4.2412*x+1.09530.4,0.5-1.5804*x3+3.4792*x2-2.7683*x+0.849840.5,0.6-0.86962*x3+2.1999*x2-2.0008*x+0.696320.6,0.7-0.46379*x3+1.3477*x2-1.4042*x+0.557120.7,0.8-0.27191*x3+0