三次样条差值

1计算原理和程序流程图

设在区间[a, b]上给定n+1个节点xi（a≤x₀<x₁< …<x_n≤b）,在节点x_i处的函数值为y_i=f(x_i) (i= 0,1,…,n)。若函数S(x)满足以下三个条件：

(1) 在每个子区间[x_i-1, x_i] (i= 0,1,…,n)上，S(x)是三次多项式；

(2) S(x_i) = y_i(i=0,1,…,n)；

(3) 在区间[a, b]上，S(x)的二阶导数S”(x)连续。

则称S(x)为函数y_i= f(x)在区间[a,b]上的三次样条插值函数。

由定义可知S(x)共有4n个待定参数，根据条件(3)可得如下3n-3个方程，再加上条件(2)的n+1 个方程，共有4n-2 个方程，因此还需要增加两个条件才可以确定出4n 个待定参数。所增加的条件称为边界条件或端点条件。常用的三种边界条件为：①已知 f”(a), f”(b)；②已知 f”(a)= f”(b)；③已知 f(x)是以 T=b-a 为周期的周期函数。

三次样条插值计算流程图如下：

 

图2求三次样条插值函数的程序框图

2.2程序使用说明

根据需要可在程序中直接修改函数f（x），插值多项式阶数和边界条件等即可。对于例题中的三弯矩方程组因为是严格对角占优的，所以选用追赶法求解。

2.3计算结果展示

例题为课本141页计算实习，分别把不同阶数的牛顿插值和三次插值多项式画图比较。结果得到的S(x）表达式如下：

y=69962787255772577/35184372088832- (54675*(x - 1)^3)/160004 - (38868241859560495*x)/17592186044416。

其中Netwon10插值误差EnMax =1.9156，三次样条插值误差EsMax =0.0220。函数图像如下：

%%% ** 文件名:spline.m% % ** 创建人:***% % ** 日 期:2016.12.14% % ** 描 述:三次样条差值运算 与牛顿插值的比较% % **   第 141 页计算实习  求三次样条插值% **   Runge函数f(x)=(1+25x^2)^(-1),(-1<= x <=1)% % **  取等距节点，构造牛顿插值多项式N5（x）和N10（x）及三次样条差值S10（x）。% **  分别将三种差值多项式与f(x)的曲线画在同一个坐标系进行比较。% **参考教材：《数值分析》李乃成，梅立泉，科学出版社   clearclcclose all%%%差值多项式N5(x)N=5;i=0:N;fprintf('插值节点xi');xx=-1+2/N*i;fprintf('插值节点xi处的函数值f(xi)');fx=(1+25*xx.^2).^(-1);yy=fx;%----------求Newton插值多项式的各阶差商-------%注意：数组下表为1...N+1，区别于课本中的0...N  在MATLAB中数组下标是从1开始的for k=1:N  %求差商表    for i=N+1:-1:k+1%之所以采用倒序，是为了不覆盖掉上一次循环所得的值        yy(i)=(yy(i)-yy(i-1))/(xx(i)-xx(i-k));    endendsyms x y%--------用秦九韶算法计算多项式的表达式y=yy(N+1);%差商表的最后一个值for i=N:-1:1    y=y*(x-xx(i))+yy(i);%之所以用逆序是为了从项数最多的到项数最少的，进而得到表达式endfprintf('Newton插值多项式的表达式:');simplify(y);%-------计算xk(k=1...99)对应的y(xk),N(xk),S(xk)------K=99;k=1:K;xk=-1+0.02*k;yk=1./(1+25*xk.^2);%-------计算N(x)在对应数据点上的的函数------Nk5=zeros(1,K);for k=1:K    Nk5(k)=yy(N+1);    for i=N:-1:1        Nk5(k)=Nk5(k)*(xk(k)-xx(i))+yy(i);%得到表达式 在k的值    endend%%%差值多项式N10(x).S10(x)N=10;i=0:N;%format short;fprintf('插值节点xi');xx=-1+2/N*ifprintf('插值节点xi处的函数值f(xi)');fx=(1+25*xx.^2).^(-1)yy=fx;%----------求Newton插值多项式的各阶差商-------%注意：数组下表为1...N+1，区别于课本中的0...N  在MATLAB中数组下标是从1开始的for k=1:N  %求差商表    for i=N+1:-1:k+1%之所以采用倒序，是为了不覆盖掉上一次循环所得的值        yy(i)=(yy(i)-yy(i-1))/(xx(i)-xx(i-k));    endendfprintf('Newton插值多项式的各阶差商');yy %最终所得的对角线上的N+1个值syms x y%--------用秦九韶算法计算多项式的表达式y=yy(N+1);%差商表的最后一个值for i=N:-1:1    y=y*(x-xx(i))+yy(i);%之所以用逆序是为了从项数最多的到项数最少的，进而得到表达式endfprintf('Newton插值多项式的表达式:');simplify(y)%----------（自然）三次样条插值多项式-----------M0=Mn=0M=zeros(1,N+1);%M(i)对应s(x)在节点x(i)处的二阶导数h=zeros(1,N);%h(i)=x(i)-x(i-1)A=zeros(N-1,N-1);%三弯矩方程组的系数矩阵（自然三次样条插值）b=zeros(1,N-1);a=b;c=b;d=b;h(1)=xx(2)-xx(1);for i=1:N-1 %i代表各数组中的第i个元素，对应的y的下标为i-1，h为i，d为i+1    h(i+1)=xx(i+2)-xx(i+1);    a(i)=h(i)/(h(i)+h(i+1));%三弯矩方程组的第-1条对角线元素  miou    c(i)=1-a(i);%三弯矩方程组的第-1条对角线元素  lanmuda    b(i)=2;%三弯矩方程组的主对角线上的元素  主对角    d(i+1)=6/(h(i)+h(i+1))*((fx(i+2)-fx(i+1))/h(i+1)-(fx(i+1)-fx(i))/h(i));endd=d(2:N);M(1)=0;M(N+1)=0;%自然三次样条插值函数的边界条件%-----------用追赶法解严格三对角占优阵----------z=zeros(N-1,1);w=z;%z为中间矩阵，w为方程的解u=zeros(1,N-1);l=u;%u矩阵U的主对角线元素，l矩阵L的第-1条对角线元素%---------追过程--------u(1)=b(1);z(1)=d(1);for i=2:N-1    l(i)=a(i)/u(i-1);    u(i)=b(i)-l(i)*c(i-1);    z(i)=d(i)-l(i)*z(i-1);end%---------赶过程--------w(N-1)=z(N-1)/u(N-1);for i=N-2:-1:1    w(i)=(z(i)-c(i)*w(i+1))/u(i);endfprintf('三弯矩方程组的解向量：');w'M(2:N)=w;for i=1:N    syms x y;    fprintf('在第%d个区间(%5.4f<=x<=%5.4f)内，S(x)的表达式为：\n',i,xx(i),xx(i...+1));    y=(xx(i+1)-x)^3/(6*h(i))*M(i)+(x-xx(i))^3/(6*h(i))*M(i+1)+(yy(i)-h(i)...^2/6*M(i))*(xx(i+1)-x)/h(i)+(yy(i+1)-h(i)^2/6*M(i+1))*(x-xx(i))/h(i);    y=simplify(y)end%%%开始准备作图%-------计算xk(k=1...99)对应的y(xk),N(xk),S(xk)------K=99;k=1:K;xk=-1+0.02*k;yk=1./(1+25*xk.^2);%-------计算N(x)在对应数据点上的的函数------Nk10=zeros(1,K);for k=1:K    Nk10(k)=yy(N+1);    for i=N:-1:1        Nk10(k)=Nk10(k)*(xk(k)-xx(i))+yy(i);%得到表达式 在k的值    endend%-------计算S(x)在对应数据点上的的函数------Sk10=zeros(1,K);s=zeros(1,K);for k=1:K    for i=1:N-1        if xk(k)<=xx(i+1)            s(k)=i;            break;      %s为xk(k)所在的插值区间号        else            s(k)=N;        end    end      hh=xx(s(k)+1)-xx(s(k));    xRight=xx(s(k)+1)-xk(k);    xLeft=xk(k)-xx(s(k));    Sk10(k)=(M(s(k))*xRight^3/6+M(s(k)+1)*xLeft^3/6+(fx(s(k))-M(s(k))*hh^2/6)...*xRight+(fx(s(k)+1)-M(s(k)+1)*hh^2/6)*xLeft)/hh;%得到表达式 在k的值endfigure(1);plot(xk,yk,'r-',xk,Nk5,'m',xk,Nk10,'g',xk,Sk10,'b');title('差值结果');legend('原函数f(x)','Newton插值函数N5(x)','Newton插值函数N10(x)','三次样条插值函数S10(x)');hold onEn=abs(Nk10-yk);fprintf('Netwon10插值误差');EnMax=max(En)Es=abs(Sk10-yk);fprintf('三次样条插值误差');EsMax=max(Es)