数值分析实验报告-插值、三次样条(教育教学)

实验报告：牛顿差值多项式&三次样条
问题：在区间[-1,1]上分别取n=10、20用两组等距节点对龙格函数2
1()25f x x 作多项式插值及三次样条插值，对每个n 值，分别画出插值函数及()f x 的图形。

实验目的：通过编程实现牛顿插值方法和三次样条方法，加深对多项式插值的理解。

应用所编程序解决实际算例。

实验要求：
1． 认真分析问题，深刻理解相关理论知识并能熟练应用；
2． 编写相关程序并进行实验；
3． 调试程序，得到最终结果；
4． 分析解释实验结果；
5． 按照要求完成实验报告。

实验原理：
详见《数值分析 第5版》第二章相关内容。

实验内容：
（1）牛顿插值多项式
1.1 当n=10时：
在Matlab 下编写代码完成计算和画图。

结果如下：
代码：
clear all
clc
x1=-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));
end
end
syms 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);
end
syms P
P=sum(p);
P10=vpa(expand(P),5);
x0=-1:0.001:1;
y0=subs(P,x,x0);
y2=subs(1/(1+25*x^2),x,x0);
plot(x0,y0,x0,y2)
grid on
xlabel('x')
ylabel('y')
P10即我们所求的牛顿插值多项式，其结果为：P10（x）=-220.94*x^10+494.91*x^8-9.5065e-14*x^7-381.43*x^6-8.504e-14*x^5+123.36*x^4+2.0202e-1 4*x^3-16.855*x^2-6.6594e-16*x+1.0
并且这里也能得到该牛顿插值多项式的在[-1，1]上的图形，并和原函数进行对比（见Fig.1）。

Fig.1 牛顿插值多项式（n=10）函数和原函数图形
从图形中我们可以明显的观察出插值函数在两端点处发生了剧烈的波动，产生了极大的误差，即龙格现象，当n=20时，这一现象将更加明显。

1.2 当n=20时：
对n=10的代码进行修改就可以得到n=20时的代码。

将“x1=-1:0.2:1;”改为“x1=-1:0.1:1;”即可。

运行程序，我们得到n=20时的牛顿插值多项式，结果为：P20(x)= 260188.0*x^20 - 1.0121e6*x^18 + 2.6193e-12*x^17 + 1.6392e6*x^16 + 2.248e-11*x^15 - 1.4429e6*x^14 - 4.6331e-11*x^13 + 757299.0*x^12 + 1.7687e-11*x^11 - 245255.0*x^10 + 2.1019e-11*x^9 + 49318.0*x^8 + 3.5903e-12*x^7 - 6119.2*x^6 - 1.5935e-12*x^5 + 470.85*x^4 + 1.3597e-14*x^3 - 24.143*x^2 - 1.738e-14*x + 1.0
同样的，这里得到了该牛顿插值多项式的在[-1，1]上的图形，并和原函数进行对比（见Fig.2）。

Fig.2牛顿插值多项式（n=20）函数和原函数图形
当n=20时，端点处发生了更加剧烈的震荡。

表明随着分段不断增加，牛顿插值多项式与原函数的误差不但没有减少，反而变得更大了。

（2）三次样条
2.1 当n=10时：
在Matlab下编写代码完成计算和画图。

结果如下：
代码：
clear all
clc
x1=-1:0.2:1;
y1=1./(1+25.*x1.^2);
syms x
m1=subs(diff(1/(1+25*x^2)),-1);
m2=subs(diff(1/(1+25*x^2)),1);
n=length(x1);
syms a b h f d
for i=1:n-1
h(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
a(i)=h(i-1)/(h(i-1)+h(i));
b(i)=h(i)/(h(i-1)+h(i));
end
d(1)=6/h(1)*(f(1)-m1);
d(n)=6/h(n-1)*(m2-f(n-1));
for i=2:n-1
d(i)=6*(f(i)-f(i-1))/(h(i-1)+h(i));
end
D=d';
A=2.*eye(n);
for i=1:n-1
A(i,i+1)=b(i);
A(i+1,i)=a(i+1);
end
M=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); end
S=vpa(expand(s.'),5);
for i=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 on
end
y2=subs(1/(1+25*x^2),x,-1:0.001:1);
plot(-1:0.001:1,y2,'r')
grid on
xlabel('x')
ylabel('y')
S即为我们所求的三次样条，其结果为：
S10(x) =
0.08225*x^3+0.36953*x^2+0.56627*x+0.31745 [-1,-0.8]
0.96279*x^3+2.4828*x^2+2.2569*x+0.76829 [-0.8,-0.6]
0.81773*x^3+2.2217*x^2+2.1002*x+0.73696 [-0.6,-0.4]
13.413*x^3+17.336*x^2+8.1461*x+1.5431 [-0.4,-0.2]
-54.471*x^3-23.394*x^2-1.8741e-17*x+1.0 [-0.2,0]
54.471*x^3-23.394*x^2+1.9683e-17*x+1.0 [0,0.2]
-13.413*x^3+17.336*x^2-8.1461*x+1.5431 [0.2,0.4]
-0.81773*x^3+2.2217*x^2-2.1002*x+0.73696 [0.4,0.6]
-0.96279*x^3+2.4828*x^2-2.2569*x+0.76829 [0.6,0.8]
-0.08225*x^3+0.36953*x^2-0.56627*x+0.31745 [0.8,1] 并且这里可以得到该三次样条的在[-1，1]上的图形，并和原函数进行对比（见Fig.3）。