利用MATLAB验证样条插值的收敛性我要分享

问题提出

多项式插值是不收敛的，即插值的节点多，效果不一定就好。对样条函数插值又如何呢？理论上证明样条插值的收敛性是比较困难的，但通过本实验可以验证这一理论结果。.

实验内容

请按一定的规则分别选择等距或者非等距的插值节点，并不断增加插值节点的个数。考虑实验2.1中的函数或选择其他你有兴趣的函数。

实验要求

（1）随节点个数增加，比较被逼近函数和样条插值函数误差的变化情况。分析所得结果并与拉格朗日多项式插值比较(可以用MATLAB的函数“spline”作此函数的三次样条插值，取n=10、20，分别画出插值函数及原函数的图形)。

（2）样条插值的思想是早产生于工业部门。作为工业应用的例子考虑如下问题：某汽车制造商用三次样条插值设计车门的曲线，其中一段的数据如下：

要求：i.自己编程计算（用三弯矩、三转角方程均可）

      ii.主函数myspline(x,y,边界类型，边界值，x_i )

        其中：x 节点   y 节点上的函数值  x_i  未知节点

        返回：S(x_i)

         iii.三对角方程组用追赶法求解(书P160)。

实验步骤及结果分析

1. 用MATLAB的函数“spline”作此函数的三次样条插值，取n=10、20，分别画出插值函数及原函数的图形。随节点个数增加，比较被逼近函数和样条插值函数误差的变化情况。分析所得结果并与拉格朗日多项式插值比较。逼近函数选用实验2.1中函数。

结果：

①原函数与三次样条插值函数图像：（程序：test2_2_1_1.m）

Elapsed time is 0.147261 seconds.

②原函数与拉格朗日插值函数图像：（程序：test2_2_1_2.m ）

Elapsed time is 5.914741 seconds.

分析：

第一幅图使用三次样条插值函数spline插值后画出的图，从图中可以看出，随着插值节点的增高，插值函数与原函数越来越逼近，其误差越来越小，而且即使在n较小时插值函数与原函数的逼近程度也很好。

第二幅图是使用拉格朗日型插值做出的插值函数与原函数的图像。从图中可以看出，插值函数与原函数在插值区域边界出现偏离，而且随着n的增大，偏离越加剧，出现了龙格现象。当n→∞时， L_n(x)不收敛于f(x),所以随着n的增大，高次拉格朗日型插值的插值函数与原函数f(x)的逼近效果不好。

此外从两种方法插值所用的时间分别为：spline插值用时t1=0.147261 s，而拉格朗日插值用时t2=5.914741 s。可以明显的看出三次样条插值要比拉格朗日插值收敛速度快许多。更重要的是，三次样条插值有更高的光滑性，所以在实际应用中三次样条插值要比拉格朗日插值应用广泛的多。

2.利用三弯矩方法编制三次样条差值主函数(myspline.m)，三对角方程组用追赶法求解(zhuigan.m)。myspline函数的输入有：x为已知节点，y为已知节点的函数值；y11、y1N分别为两端点的一阶导数值，x0为未知节点；返回未知节点的函数值y0及所在区间的三次样条插值多项式。

结果：

运行文件test2_2_2.m后就可以得到插值后的数据为：

附录

附录一：myspline.m文件源码：

%用三弯矩方法编制的三次样条插值函数
%  x为已知节点
%  y为已知节点的函数值
%  y11、y1N分别为两端点的一阶导数值，
%  x0为未知节点；
%  返回未知节点的函数值y0及所在区间的三次样条插值多项式。

function [S,y0]=myspline(x,y,y11,y1N,x0)
format short
syms w;
n=length(x)-1;
A=diag(2*ones(1,n+1));
u=ones(n,1);
lamda=ones(n,1);
d=zeros(n+1,1);
for i=1:n-1
    u(i)=(x(i+1)-x(i))/(x(i+2)-x(i));
    lamda(i+1)=(x(i+2)-x(i+1))/(x(i+2)-x(i));
    d(i+1)=6*(((y(i+2)-y(i+1))*(x(i+1)-x(i))-(y(i+1)-y(i))*(x(i+2)-x(i+1)))/((x(i+2)-x(i+1))*(x(i+1)-x(i))*(x(i+2)-x(i))));
end
d(1)=6/(x(2)-x(1))*((y(2)-y(1))/(x(2)-x(1))-y11);
d(n+1)=6/(x(n+1)-x(n))*(y1N-(y(n+1)-y(n))/(x(n+1)-x(n)));
for i=1:n
    A(i+1,i)=u(i);
    A(i,i+1)=lamda(i);
end
M=zhuigan(A,d);
for i=1:n
    if(x(i)=x0)
        j=i;
        break;
    end
end
h=x(j+1)-x(j);
S=M(j)*((x(j+1)-w)^3/(6*h))+M(j+1)*((w-x(j))^3/(6*h))+(y(j)-M(j)*h*h/6)*(x(j+1)-w)/h+(y(j+1)-M(j+1)*h*h/6)*(w-x(j))/h;
y0=subs(S,w,x0);

附录二：test2_2_1_1.m文件源码：

%本次使用的是三次样条差值。
tic
for n=10:10:20;    %插值次数n=10、20
x=-1:2/n:1;        %等距节点
y=1./(1+25*x.^2);  %原函数   
x1=-1:0.01:1;      %画图取得离散点
y1=spline(x,y,x1); %三次样条差值函数  直接使用spline函数
f=1./(1+25*x1.^2); %原函数,'linewidth',4,'linewidth',2
subplot(1,2,n/10)
plot(x1,f,'b-');hold on;
plot(x1,y1,'r-');hold off;
title(['原函数与',num2str(n),'次三次样条差值函数图像']);
xlabel('x');ylabel('y');
legend('原函数','插值函数');
grid on
end
toc

附录三：test2_2_1_2.m文件源码： 

%本次使用的是三次样条差值。
tic
for n=10:10:20;            %n为插值阶数
subplot(1,2,n/10)          %将各阶原函数与插值函数的图像分别一张图片上面。
syms x;
f=1/(1+25*x^2);            %原函数  
x1=sym(zeros(n+1));
W=sym(ones(n+1));
L=sym(0);
for i=0:n
    x1(i+1)=-1+2*i/n;      %划分等距节点
end
for j=0:n
    for i=0:n
     if i~=j               %求插值基函数
         w=(x-x1(i+1))/(x1(j+1)-x1(i+1));   
         W(j+1)=W(j+1)*w;
     end
    end
      L=L+W(j+1)*(1/(1+25*x1(j+1)^2));     %插值函数
end
LL(n)=simplify(L);        %插值函数化简
x=-1:0.01:1;
y1=subs(f,x);             %将x带入原函数f中
y2=subs(L,x);             %将x带入插值函数L中
plot(x,y1,'b'); hold on;
plot(x,y2,'r'); hold off;
title(['原函数与',num2str(n),'次插值函数']);
xlabel('x'); ylabel('y');
legend('原函数','插值函数');
grid on
end
toc

附录四：test2_2_2.m文件源码： 

 x=0:10;
 y=[0.0 0.79 1.53 2.19 2.71 3.03 3.27 2.89 3.06 3.19 3.29];
 for x0=0:10
    [S,y0]=myspline(x,y,0.8,0.2,x0);
    S1=diff(S);
    y1(x0+1)=subs(S1,w,x0);
 end
y1

附录五：zhuigan.m文件源码： 

function x=zhuigan(A,f)
n=rank(A);
b=ones(n,1);
a=ones(n-1,1);
c=ones(n-1);
for i=1:n-1
    a(i,1)=A(i+1,i);
    c(i,1)=A(i,i+1);
    b(i,1)=A(i,i);
end
b(n,1)=A(n,n);
d(1,1)=c(1,1)/b(1,1);
for i=2:n-1
    d(i,1)=c(i,1)/(b(i,1)-a(i-1,1)*d(i-1,1));
end
y(1,1)=f(1,1)/b(1,1);
for i=2:n
    y(i,1)=(f(i,1)-a(i-1,1)*y(i-1,1))/(b(i,1)-a(i-1,1)*d(i-1,1));
end
x(n,1)=y(n,1);
for i=(n-1):-1:1
    x(i,1)=y(i,1)-d(i,1)*x(i+1,1);
end

注意：以上几个m文件放在同一目录下。