三次样条插值(附上MATLAB代码)

三次样条插值(附上MATLAB代码)

前言：

工程数值方法课程上，老师讲了拉格朗日插值法、哈密顿插值法以及三次样条插值，但是三次样条插值实在难以理解，故写一篇博客，以后理解也方便。

问题起始：

龙格现象

使用均匀节点构造高次多项式时，由于节点数量不够，导致在边缘会出现误差极大的情况。

此时，有一种优化思路便为不考虑多个节点，而考虑相邻两个节点的影响，利用相邻两个节点的相关性质来寻找一条多项式进行插值。这里选择使用三次多项式（一次二次不够拟合，三次以上即可，高次多项式可以使用相似的条件。

基本思路：

核心思路：

对于给定点(x0,y0),……，(xn,yn),在每一个子区间[xi,xi+1]中，局部构建出一个三次曲线来进行拟合。此时，三次曲线的形式可以写成如下形式：
$$S_i\left(x\right)=a_i+b_i\left(x-x_i\right)+c_i{\left(x-x_{\mathrm{i}}\right)}^2+{\mathrm{d}}_{\mathrm{i}}{\left(\mathrm{x}-{\mathrm{x}}_{\mathrm{i}}\right)}^3$$
解这个三次多项式，至少需要四个条件。它满足如下四个条件：
$$\{\begin{array}{l}{\mathrm{S}}_{\mathrm{i}}\left({\mathrm{x}}_{\mathrm{i}}\right)={\mathrm{y}}_{\mathrm{i}},{\mathrm{S}}_{\mathrm{i}}\left({\mathrm{x}}_{\mathrm{i}+1}\right)={\mathrm{y}}_{\mathrm{i}+1}\\ \mathrm{S}{\prime }_{\mathrm{i}-1}\left({\mathrm{x}}_{\mathrm{i}}\right)=\mathrm{S}{\prime }_{\mathrm{i}}\left({\mathrm{x}}_{\mathrm{i}}\right)\\ {\mathrm{S}}_{\mathrm{i}-1}^{\prime \prime }\left({\mathrm{x}}_{\mathrm{i}}\right)={\mathrm{S}}_{\mathrm{i}}^{\prime \prime }\left({\mathrm{x}}_{\mathrm{i}}\right)\end{array}$$
这四个条件的含义为，对于该曲线，两端由两点确定，对于第i个曲线，满足前后的两条曲线在xi点的变化程度相等。但是注意，这里埋了一个小坑，在端点处，两端的斜率的条件是无法计算的。所以实际求解时，需要对端点处进行条件补充。

求解关系：

此时，根据这三个限制可以进行求解。我们假定在[xi,xi+1]这个区间上进行求解，则有：

(1) 约束一：

ai此时由第一个条件，可以直接用yi来代替，利用第二个条件可以得到：
$${\mathrm{S}}_{\mathrm{i}}\left({\mathrm{x}}_{\mathrm{i}+1}\right)={\mathrm{y}}_{\mathrm{i}+1}$$

$${\mathrm{y}}_{\mathrm{i}+1}={\mathrm{y}}_{\mathrm{i}}+{\mathrm{b}}_{\mathrm{i}}\left(\mathrm{x}-{\mathrm{x}}_{\mathrm{i}}\right)+{\mathrm{c}}_{\mathrm{i}}{\left(\mathrm{x}-{\mathrm{x}}_{\mathrm{i}}\right)}^2+{\mathrm{d}}_{\mathrm{i}}{\left(\mathrm{x}-{\mathrm{x}}_{\mathrm{i}}\right)}^3$$

为方便起见，我们令：
$$\{\begin{array}{l}{\Delta }_{\mathrm{i}}={\mathrm{y}}_{\mathrm{i}+1}-{\mathrm{y}}_{\mathrm{i}}\\ {\delta }_{\mathrm{i}}={\mathrm{x}}_{\mathrm{i}+1}-{\mathrm{x}}_{\mathrm{i}}\end{array}$$
此时则有式子一：
$$\frac{{\Delta }_{\mathrm{i}}}{{\delta }_{\mathrm{i}}}={\mathrm{b}}_{\mathrm{i}}+{\delta }_{\mathrm{i}}{\mathrm{c}}_{\mathrm{i}}+{{\delta }_{\mathrm{i}}}^2{\mathrm{d}}_{\mathrm{i}}$$
(2)约束二：
$$\mathrm{S}{\prime }_{\mathrm{i}-1}\left({\mathrm{x}}_{\mathrm{i}}\right)=\mathrm{S}{\prime }_{\mathrm{i}}\left({\mathrm{x}}_{\mathrm{i}}\right)$$
由于式子中具有x-xi项，对于右边求导后代入不难得到：
$$\mathrm{S}{\prime }_{\mathrm{i}}\left({\mathrm{x}}_{\mathrm{i}}\right)={\mathrm{b}}_{\mathrm{i}}$$
此时左边为：
$$\begin{gathered} \mathrm{S}{\prime }_{\mathrm{i}-1}\left(\mathrm{x}\right)={\mathrm{b}}_{\mathrm{i}-1}+2{\mathrm{c}}_{\mathrm{i}-1}\left(\mathrm{x}-{\mathrm{x}}_{\mathrm{i}-1}\right)+3{\mathrm{d}}_{\mathrm{i}-1}{\left(\mathrm{x}-{\mathrm{x}}_{\mathrm{i}-1}\right)}^2 \\ \mathrm{S}{\prime }_{\mathrm{i}-1}\left({\mathrm{x}}_{\mathrm{i}}\right)={\mathrm{b}}_{\mathrm{i}-1}+2{\mathrm{c}}_{\mathrm{i}-1}{\delta }_{\mathrm{i}-1}+3{\mathrm{d}}_{\mathrm{i}-1}{{\delta }_{\mathrm{i}-1}}^2=b_i \end{gathered}$$
(3)约束三思路与约束二相同，直接写出最后结果：
$${\mathrm{S}}_{\mathrm{i}-1}^{\prime \prime }\left({\mathrm{x}}_{\mathrm{i}}\right)=2{\mathrm{c}}_{\mathrm{i}-1}+6{\delta }_{\mathrm{i}-1}{\mathrm{d}}_{\mathrm{i}-1}={\mathrm{S}}_{\mathrm{i}}^{\prime \prime }\left({\mathrm{x}}_{\mathrm{i}}\right)=2{\mathrm{c}}_{\mathrm{i}}$$
此时，我们有了三个未知数以及三个约束方程，可以进行求解，求解之前先对方程进行化简，求解出其中一个未知数，另外两个未知数用这个未知数进行比较即可。考虑到将ci进行取代较为复杂，因此选择ci为求解的未知数，用ci表示di和bi，可以在约束二，约束三中进行取代，则有：
$$\{\begin{array}{l}{\mathrm{d}}_{\mathrm{i}-1}=\frac{1}{3{\delta }_{\mathrm{i}-1}}\left({\mathrm{c}}_{\mathrm{i}}-{\mathrm{c}}_{\mathrm{i}-1}\right)\to {\mathrm{d}}_{\mathrm{i}}=\frac{1}{3{\delta }_{\mathrm{i}}}\left({\mathrm{c}}_{\mathrm{i}+1}-{\mathrm{c}}_{\mathrm{i}}\right)\\ {\mathrm{b}}_{\mathrm{i}}=\frac{{\Delta }_{\mathrm{i}}}{{\delta }_{\mathrm{i}}}-\frac{2}{3}{\delta }_{\mathrm{i}}{\mathrm{c}}_{\mathrm{i}}-\frac{{\delta }_{\mathrm{i}}}{3}{\mathrm{c}}_{\mathrm{i}+1}\end{array}$$
代入约束一的方程中，整理后可写成：
$${\delta }_{\mathrm{i}-1}{\mathrm{c}}_{\mathrm{i}-1}+2\left({\delta }_{\mathrm{i}-1}+{\delta }_{\mathrm{i}}\right){\mathrm{c}}_{\mathrm{i}}+{\delta }_{\mathrm{i}}{\mathrm{c}}_{\mathrm{i}+1}=3\left(\frac{{\Delta }_{\mathrm{i}}}{{\delta }_{\mathrm{i}}}-\frac{{\Delta }_{\mathrm{i}-1}}{{\delta }_{\mathrm{i}-1}}\right)$$
注意，此时ci为未知数，此时为特殊的一种方程形式，最后整理可以化成三对角方程，三对角方程可以利用追赶法进行快速求解。追赶法的函数如下，具体原理不再赘述：

function x=tridecomposition(up,down,angle,b)
    n=length(angle);
    for i = 2:1:n
        down(i-1) = down(i-1)/angle(i-1);
        angle(i) = angle(i) - up(i-1) * down(i-1);
        b(i) = b(i) - b(i-1) * down(i-1);
    end
    x(n)=b(n)/angle(n);
    for i=n-1:-1:1
        x(i)=(b(i)-up(i)*x(i+1))/angle(i);
    end
end

此时的方程形式为：
$$\left[\begin{array}{lccccc}{\delta }_0& 2\left({\delta }_0+{\delta }_1\right)& {\delta }_1& & & \\ & {\delta }_1& 2\left({\delta }_1+{\delta }_2\right)& {\delta }_2& & \\ & & \ddots & \ddots & \ddots & \\ & & & {\delta }_{\mathrm{n}-2}& 2\left({\delta }_{\mathrm{n}-2}+{\delta }_{\mathrm{n}-1}\right)& {\delta }_{\mathrm{n}-1}\end{array}\right]\left[\begin{array}{c}{\mathrm{c}}_{\mathrm{o}}\\ {\mathrm{c}}_1\\ ⋮\\ {\mathrm{c}}_{\mathrm{n}}\end{array}\right]=3\left[\begin{array}{c}\frac{{\Delta }_1}{{\delta }_1}-\frac{{\Delta }_0}{{\delta }_0}\\ \frac{{\Delta }_2}{{\delta }_2}-\frac{{\Delta }_1}{{\delta }_1}\\ ⋮\\ \frac{{\Delta }_{\mathrm{n}-1}}{{\delta }_{\mathrm{n}-1}}-\frac{{\Delta }_{\mathrm{n}-2}}{{\delta }_{\mathrm{n}-2}}\end{array}\right]$$
此时方程左边为n-1个条件，有n+1个未知数，仍然缺两个条件，即上文提到的端点条件

端点条件：

端点条件，老师介绍了三种思路：

1. 自然样条（Natural Spline)：

$${\mathrm{S}}_0^{\prime \prime }\left({\mathrm{x}}_0\right)=0;{\mathrm{S}}_{\mathrm{n}-1}^{\prime \prime }\left({\mathrm{x}}_{\mathrm{n}}\right)=0$$

​ 评价为简单粗暴的直接赋值

2. 夹紧的三次样条(Clamped Cubic Spline)
   $$S{\prime }_0\left(x_0\right)=v_0;S{\prime }_{n-1}\left(x_n\right)=v_n$$
   需要引入额外的条件进行赋值，对自然样条的一些补充，但基本思路一致。

3. 非结三次样条（Not-a-knot cubic Spline)-MATLAB默认的spline函数命令

寻找三次导数的值进行代入，其实就是带入
$${\mathrm{S}}_0^{\prime \prime \prime }\left({\mathrm{x}}_1\right)=6{\mathrm{d}}_1;{\mathrm{S}}_{\mathrm{n}-2}^{\prime \prime \prime }\left({\mathrm{x}}_{\mathrm{n}-1}\right)=6{\mathrm{d}}_{\mathrm{n}-1}$$
这里只编写前两个（懒，作业只要求前两个qaq）

那么前两种可以转化为：

自然样条：
$$\{\begin{array}{l}{\mathrm{c}}_0=0\\ {\mathrm{c}}_{\mathrm{n}}=0\end{array}$$
加紧的三次样条：
$$\{\begin{array}{l}2{\delta }_0{\mathrm{c}}_0+{\delta }_0{\mathrm{c}}_1=3\left(\frac{{\Delta }_0}{{\delta }_0}-{\mathrm{v}}_0\right)\\ {\delta }_{\mathrm{n}-1}{\mathrm{c}}_{\mathrm{n}-1}+2{\delta }_{\mathrm{n}-1}{\mathrm{c}}_{\mathrm{n}}=3\left({\mathrm{v}}_{\mathrm{n}}-\frac{{\Delta }_{\mathrm{n}-1}}{{\delta }_{\mathrm{n}-1}}\right)\end{array}$$
补充进方程即可求解。求解出cn后就是简单的回代，然后写成多项式即可。

最终求解方程的形式：
$$\left[\begin{array}{lccccc}2{\delta }_0& {\delta }_0& & & & \\ {\delta }_0& 2\left({\delta }_0+{\delta }_1\right)& {\delta }_1& & & \\ & {\delta }_1& 2\left({\delta }_1+{\delta }_2\right)& {\delta }_2& & \\ & & \ddots & \ddots & \ddots & \\ & & & {\delta }_{\mathrm{n}-2}& 2\left({\delta }_{\mathrm{n}-2}+{\delta }_{\mathrm{n}-1}\right)& {\delta }_{\mathrm{n}-1}\\ & & & & {\delta }_{\mathrm{n}-1}& 2{\delta }_{\mathrm{n}-1}\end{array}\right]\left[\begin{array}{c}{\mathrm{c}}_{\mathrm{o}}\\ {\mathrm{c}}_1\\ ⋮\\ {\mathrm{c}}_{\mathrm{n}}\end{array}\right]=3\left[\begin{array}{c}\frac{{\Delta }_0}{{\delta }_0}-{\mathrm{V}}_0\\ \frac{{\Delta }_1}{{\delta }_1}-\frac{{\Delta }_0}{{\delta }_0}\\ \frac{{\Delta }_2}{{\delta }_2}-\frac{{\Delta }_1}{{\delta }_1}\\ ⋮\\ \frac{{\Delta }_{\mathrm{n}-1}}{{\delta }_{\mathrm{n}-1}}-\frac{{\Delta }_{\mathrm{n}-2}}{{\delta }_{\mathrm{n}-2}}\\ {\mathrm{v}}_{\mathrm{n}}-\frac{{\Delta }_{\mathrm{n}-1}}{{\delta }_{\mathrm{n}-1}}\end{array}\right]$$
标准三对角形式，函数一套结束。Perfect！！！！！！！！！！！

MATLAB代码:

调用指南：

当flag=0时，端点条件取为Natural Spline；

当flag=1时，端点条件取为Clamped Cubic Spline，若斜率不输入，默认为0）；

当flag=2时，调用Matlab自带的spline命令。

function v=myCubicSpline(x,y,u,flag,varargin)
     if(nargin==4)
        v0=0;
        vn=0;
     end
     if(nargin==5)
        v0=varargin{1};
        vn=0;
     end
     if(nargin==6)
        v0=varargin{1};
        vn=varargin{2};
     end
     %初始化
     delta=zeros(length(x)-1,1);
     Delta=zeros(length(x)-1,1);
     for i=1:1:length(x)-1
        delta(i)=x(i+1)-x(i);
        Delta(i)=y(i+1)-y(i);
     end
     angle=zeros(length(x),1);
     rightb=zeros(length(x),1);
     for i=1:1:length(x)
        if(i==1) 
            angle(i)=2*delta(i);
            rightb(i)=3*(Delta(i)/delta(i)-v0);
        end
        if(i>1)&&(i<length(x))
            angle(i)=2*(delta(i)+delta(i-1));
            rightb(i)=3*(Delta(i)/delta(i)-Delta(i-1)/delta(i-1));
        end
        if(i==length(x))
            angle(i)=2*delta(i-1);
            rightb(i)=3*(vn-Delta(i-1)/delta(i-1));
        end
     end
     if(flag==1)
     %三对角求解
        c=tridecomposition(delta,delta,angle,rightb);%含c0;
        d=zeros(length(x)-1,1);
        b=zeros(length(x)-1,1);
        for i=1:1:length(x)-1
            d(i)=(c(i+1)-c(i))/(3*delta(i));
            b(i)=Delta(i)/delta(i)-2/3*delta(i)*c(i)-delta(i)*c(i+1)/3;
        end
        v=zeros(length(u),1);
        %回代
        for i=1:1:length(u)
            for j=1:1:length(x)-1
                if(u(i)>=x(j))&&(u(i)<=x(j+1))
                    v(i)=y(j)+b(j)*(u(i)-x(j))+c(j)*(u(i)-x(j))^2+d(j)*(u(i)-x(j))^3;
                    break;
                end
            end
        end
        return;
     end
     if(flag==0)
         rightb(1)=0;
         rightb(length(x))=0;
         angle(1)=1;
         angle(length(x))=1;
         up=delta;
         up(length(x)-1)=0;
        c=tridecomposition(up,up,angle,rightb);%含c0;
        d=zeros(length(x)-1,1);
        b=zeros(length(x)-1,1);
        for i=1:1:length(x)-1
            d(i)=(c(i+1)-c(i))/(3*delta(i));
            b(i)=Delta(i)/delta(i)-2/3*delta(i)*c(i)-delta(i)*c(i+1)/3;
        end
        v=zeros(length(u),1);
        %回代
        for i=1:1:length(u)
            for j=1:1:length(x)-1
                if(u(i)>=x(j))&&(u(i)<=x(j+1))
                    v(i)=y(j)+b(j)*(u(i)-x(j))+c(j)*(u(i)-x(j))^2+d(j)*(u(i)-x(j))^3;
                    break;
                end
            end
        end
        return;
     end
     if(flag==2)
        v=spline(x,y,u);
     end
end

function x=tridecomposition(up,down,angle,b)
    n=length(angle);
    for i = 2:1:n
        down(i-1) = down(i-1)/angle(i-1);
        angle(i) = angle(i) - up(i-1) * down(i-1);
        b(i) = b(i) - b(i-1) * down(i-1);
    end
    x(n)=b(n)/angle(n);
    for i=n-1:-1:1
        x(i)=(b(i)-up(i)*x(i+1))/angle(i);
    end
end

比较图像