关于三次样条插值，计算方法比较复杂，但是静下心来仔细研究也是可以理解的。

本文借鉴文章来源：http://www.cnki.com.cn/Article/CJFDTotal-BGZD200611035.htm

定义：

     简单来说就是给定了一些在区间[a,b]的数据点{x1,x2,x3.....xn}，对应函数值{y1,y2,y3.....yn}，函数在[xj,xj+1]  (j=1,2,...n-1此处根据你的编译器所定，matlab数组下标从1开始的)上有表达式S(x)，且满足下面条件：

1. S(x)是一个三次多项式，在这里设为

2. S(xj)=yj                    （2-1）

3. S(xj-0)=S(xj+0)  (j=2,3....,n-1)这就是保证了在非端点处的其它点连续               (2-2)

4. S'(xj-0)=S'(xj+0) (j=2,3....,n-1)一阶导数连续                              (3-1)

5. S''(xj-0)=S''(xj+0) (j=2,3....,n-1)二阶导数连续                              (3-2)

【注】3 4 5的区间从2开始到n-1是因为两个端点不需要判断是否连续，端点处没连续这一说。

       有一个说法“n个未知数需要n个方程才能确定唯一解”，具体对不对，可以参考线性代数的知识。我们的最终目标是求出每个区间的式(1)或者函数值。 共有n-1个区间，每个区间四个参数aj，bj，cj，dj，那么就总共需要4(n-1)个求未知数。在(2-1)中给出了n个方程，(2-2)(3-1)(3-2)总共给出了3(n-2)个方程，所以依据唯一解方法可知还需要4(n-1)-3(n-2)-n=2个方程。

       对于剩下的两个方程，三次样条插值给出了两个边界约束方程，刚好凑齐两个，并且有三种，可以依照自己的兴趣选择一种便于实现的。

1. 给定了端点处的一阶导数值：S'(x1)=y1'    S'(xn)=yn'

2.给定了端点处的二阶导数值：S''(x1)=y1''    S''(xn)=yn''。自然边界条件：y1''=yn''=0

3.给定了一个周期性条件：如果f(x)是以b-a为周期的函数，在端点处便满足：S'(x1+0)=S'(xn-0)，S''(x1+0)=S''(xn-0)

下面的推导是以边界方程1为例的：

推导过程：

(一) 利用二阶导得到一些式子：

      S(x)为每个区间的三次多项式，那么S''(x)就是一次多项式。假设S''(xj)=Mj，S''(xj+1)=Mj+1的值已知，那么：

       

其中hj=xj+1-xj，然后利用S''(x)求S(x)：



 

（二）依据S(x)得到一些式子

按照上面的证明可以得到：

其中：hj=xj+1-xj

依据S(xj)=yj和S(xj+1)=yj+1得到：



——————更新日志2020-6-13————————

感谢评论区老哥@meng_zhi_xiang ，公式(5)的A下标不是j+1,应该是j

解得(求解过程就不写了，简单的二元方程组)

                                     （6-1）

          （6-2）

将Aj和Bj带入S(x)中可以得到S(x)的最终式子：



原文此处应该是有错误，Aj和Bj都没有xj的式子，但是原文的结果中包含。

 

这里就得到了S(x)的雏形了，xj、xj+1、hj都是已知的，但是Mj和Mj+1是假设已知的，下面就是求它们了。

（三）利用一阶导数得到一些式子

依据式(7)求出一阶导：

然后为了使在xj处连续平滑，那么在xj处的一阶导必须相等。即要满足S'(xj-0)=s'(xj+0)

即等式坐标表示S(x)在[xj-1,xj]的xj的一阶导值，右边表示[xj,xj+1]的xj一阶导值

其中利用S(x)的表达式可以得到等式左右两边的值：





 

由上得到两个式子：





 

由两式相等整理可得：



令：

则μjMj-1+2Mj+(1-μj)Mj+1=dj（j=2,3,...,n-1）

(四)带入边界条件

  此处选择边界条件1，即

计算：

结果写为2M1+M2=β1

 

结果写为Mn-1+2Mn=βn

（五）结果

结合（三）、（四）可以看出所有的n个式子已经齐全：

包括（三）结尾算得μjMj-1+2Mj+(1-μj)Mj+1=dj（j=2,3,...,n-1）的n-2个式子加上(四)得到的两个式子，刚好n个式子.

用矩阵表示出这n个式子即为：



方程中

   

并且hj=xj+1-xj

这样便解出了矩阵M，然后带入S(x)的式子即上面算得：



这样便求得了每一个区间上的S(x)了。

 

matlab代码：

SplineThree.m



<pre name="code" class="plain">function f = SplineThree(x,y,dy1,dyn,x0)
%x，y为输入的已知点
%x0为待求插值点的横坐标
%dy1,dyn为约束条件，是端点处的一阶导数值
n=length(x);
for j=1:n-1
    h(j)=x(j+1)-x(j);
end

%得到式子右边的结果矩阵
d(1,1)=6*((y(2)-y(1))/h(1)-dy1)/h(1);   %等式(11)的结果矩阵的上端点值
d(n,1)=6*(dyn-(y(n)-y(n-1))/h(n-1))/h(n-1); %等式(11)的结果矩阵的下端点值
for i=2:n-1
    u(i)=h(i-1)/(h(i-1)+h(i));
    d(i,1)=6*((y(i+1)-y(i))/h(i)-(y(i)-y(i-1))/h(i-1))/(h(i-1)+h(i));
end

%得到系数矩阵
A(1,1)=2;
A(1,2)=1;
A(n,n-1)=1;
A(n,n)=2;
for i=2:n-1
    A(i,i-1)=u(i);
    A(i,i)=2;
    A(i,i+1)=1-u(i);
end

%解方程
M=A\d;

format long
syms t;
%得到每个区间的式子
for i=1:n-1
   a(i)=y(i+1)-M(i+1)*h(i)^2/6-((y(i+1)-y(i))/h(i)-(M(i+1)-M(i))*h(i)/6)*x(i+1);
   b(i)=((y(i+1)-y(i))/h(i)-(M(i+1)-M(i))*h(i)/6)*t;
   c(i)=(t-x(i))^3*M(i+1)/(6*h(i));
   e(i)=(x(i+1)-t)^3*M(i)/(6*h(i));
   f(i)=a(i)+b(i)+c(i)+e(i);
    %f(i)=M(i)*(x(i+1)-t)^3/(6*h(i))+M(i+1)*(t-x(i))^3/(6*h(i))+(y(i)-M(i)*h(i)^2/6)*(x(i+1)-t)/h(i)+(y(i+1)-x(i+1)*h(i)^2/6)*(t-x(i))/h(i);
    % f(i)=((x(j+1)-x)^3)*M(i)/(6*h(i))+((x-x(i))^3)*M(i+1)/(6*h(i))+(y(i)-M(i)*(h(i)^2)/6)*((x(i+1)-x)/h(i))+(y(i+1)-(M(i+1)*(h(i)^2)/6))*((x-x(i))/h(i));
end

%化简
 f=collect(f);
 f=vpa(f,6);

 
if(nargin==5)
   nn=length(x0);
for i=1:nn
    for j=1:n-1
        if(x0(i)>=x(j)&x0(i)<=x(j+1))
             yynum(i)=subs(f(j),'t',x0(i));   %计算插值点的函数值.subs是替换函数，把x0用t替换
        end
    end
end   
f=yynum;
else
    f=collect(f);          %将插值多项式展开
    f=vpa(f,6);            %将插值多项式的系数化成6位精度的小数
end
end


 


 


SplineThree.m



<pre name="code" class="plain">% x=[-1.5 0 1 2];
% y=[0.125 -1 1 9];
% dy1=0.75;
% dyn=14;
% x0=1.5;
% answer=SplineThree(x,y,dy1,dyn);
% 
% X=[-1.5 0 1 2];
% Y=[0.125 -1 1 9];
% dY=[0.75 14]
% m=5;
% spline3(X,Y,dY,x0,m)

X=0:2*pi;
Y=sin(X);
dY=[1,1];
dy1=1;
dyn=1;
xx=0:0.5:6;
m=5;
% for i=1:length(xx)
% spline3(X,Y,dY,xx,m);
% end
yy=SplineThree(X,Y,dy1,dyn,xx);
plot(xx,yy,'r')

本文已经同步到微信公众号中，公众号与本博客将持续同步更新运动捕捉、机器学习、深度学习、计算机视觉算法，敬请关注



 

	

内容提要：三次样条插值的定义、三次样条插值函数的建立方法、不同边值条件下求解三次样条插值函数的方法.
《递推算法与多元插值》简介
本书详细介绍了多元差商与多元逆差商的递推算法及其在多元多项式与连分式插值中的应用。内容包括常用的张量积型二元多项式与连分式插值方法概述、直角三点组上的二元多项式与连分式插值及其比较研究、直角三点组上二元多项式插值余项等的进一步研究、非矩形网格上的二元多项式插值、基于二元递推多项式的散乱数据插值、基于二元连分式的散乱数据插值递推格式、非张量积型二元连分式插值、金字塔型网格点上的三元分叉连分式插值等。
本书可作为计算数学、应用数学等学科高年级本科生、硕士生、博士生数值分析、数值逼近、计算几何等相关课程的参考书或专业著作，还可供从事数值逼近与计算几何、计算机辅助几何设计、图像处理及相关领域的科技工作者参考。
本书获国家科学技术学术著作出版基金资助，入选“信息与计算科学丛书”。

所谓三次样条插值对于一个区间(a,b)将区间分成x0 = a < x1 ......xn-1 < b = xn 的n-1个区间，我们需要通过已知的n+1个点来模拟一个未知的函数，在三次样条插值中我们采用分段的方法来做这件事情。
三次样条插值得到的分段函数保证一下条件成立，而这些条件也是用来求解每一段样条插值的条件：
1 模拟出来的函数在已知点的函数值等于f的函数值
2模拟出来的分段函数是二阶连续的也就是说导数和二阶导数在分段的交界点是相等的
3需要知道在a和b点的二阶导数的情况，或者二阶导数在这n+1个点的变化规律
下面直接转载其内在的规律
已知:
a. n+1个数据点[xi, yi], i = 0, 1, …, n
b. 每一分段都是三次多项式函数曲线
c. 节点达到二阶连续
d. 左右两端点处特性(自然边界，固定边界，非节点边界)
根据定点，求出每段样条曲线方程中的系数，即可得到每段曲线的具体表达式。
插值和连续性：
, 其中 i = 0, 1, …, n-1
微分连续性：
 , 其中 i = 0, 1, …, n-2
样条曲线的微分式：
将步长
 带入样条曲线的条件：
a. 由 (i = 0, 1, …, n-1)推出
b. 由 (i = 0, 1, …, n-1)推出
c. 由  (i = 0, 1, …, n-2)推出
由此可得：
d. 由  (i = 0, 1, …, n-2)推出
设 ，则
a.  可写为：
 ，推出
b. 将ci, di带入  可得：
c. 将bi, ci, di带入 (i = 0, 1, …, n-2)可得：
这样我们可以构造一个以m为未知数的线性方程组，而且在端点条件已知的情况下我们是知道其中几个未知数的值的
端点条件
由i的取值范围可知，共有n-1个公式， 但却有n+1个未知量m 。要想求解该方程组，还需另外两个式子。所以需要对两端点x0和xn的微分加些限制。 选择不是唯一的，3种比较常用的限制如下。
a. 自由边界(Natural)
首尾两端没有受到任何让它们弯曲的力，即 。具体表示为 和 
则要求解的方程组可写为：
b. 固定边界(Clamped)
首尾两端点的微分值是被指定的，这里分别定为A和B。则可以推出
将上述两个公式带入方程组，新的方程组左侧为
c. 非节点边界(Not-A-Knot)
指定样条曲线的三次微分匹配，即
根据 和 ，则上述条件变为
新的方程组系数矩阵可写为：
右下图可以看出不同的端点边界对样条曲线的影响：
1.3 算法总结
假定有n+1个数据节点
a. 计算步长 (i = 0, 1, …, n-1)
b. 将数据节点和指定的首位端点条件带入矩阵方程
c. 解矩阵方程，求得二次微分值。该矩阵为三对角矩阵，具体求法参见我的上篇文章：三对角矩阵的求解。
d. 计算样条曲线的系数：
其中i = 0, 1, …, n-1
e. 在每个子区间 中，创建方程
原创作者博客地址：https://blog.csdn.net/flyingleo1981/article/details/53008931

	

三次样条(cubic spline)插值
1. 拟合与插值2. 分段三次样条插值3. 三对角线性系统的求解4. 固定边界三次样条的求解5. 附录5.1 龙格库塔现象
1. 拟合与插值
一个情形
有时我们获得了某条曲线上若干个离散点的信息，却不知道具体的曲线方程，又希望得到这些离散点之外的信息。
两种方法
拟合：不要求方程通过所有已知点，讲究神似，曲线整体趋势一致即可。
插值：要求方程通过所有已知点，讲究形似。
样条插值
在数值分析这个数学分支中，样条插值是使用一种名为样条的特殊分段多项式进行插值的形式。由于样条插值可以使用低阶多项式样条实现较小的插值误差，这样就避免了使用高阶多项式所出现的龙格现象(见附录1)，所以样条插值得到了流行。常用的是三次样条插值。
2. 分段三次样条插值
所谓分段三次样条插值，顾名思义，就是将区间  分成  个区间  ，共有  个点，其中端点值  。三次样条就是说每个分段区间上的曲线都是一个三次多项式方程  ，并满足以下条件。
1. 在每个分段区间  上， 都是一个三次多项式方程。
2. 在所有已知点上满足插值条件，即  。
3. 曲线光滑，即  平滑， 连续。
我们可将构造如下的三次多项式方程：
我们称这个方程为三次样条函数  。显然，每个分段区间上的三次样条函数均有  个未知系数，有  个分段区间，共有  个未知系数，需要  个方程来求解。
分析
我们需要寻找出  个方程去求解  个未知系数。
1. 所有已知点满足插值条件，即  。除端点外的  内部点  满足  ，共有  个方程；两个端点满足  ，共有  个方程。总计  个方程。
2.  个内部点一阶导数连续，即  。总计  个方程。
3.  个内部点二阶导数连续，即  。总计  个方程。
这样，我们就找出了  个方程了，还差  个方程就可以解出所有未知系数了。这  个方程我们通过边界条件给出。
自然边界(natural spline)：指定端点的二阶导数为  ，即  。
固定边界(clamped spline)：指定端点一阶导数为特定值，即  。
非扭结边界(not-a-knot spline)：使最前面两个插值点的三阶导数值相等、最后面两个插值点的三阶导数值相等，即  。
推导
为了方便推导，我们令：
由  可得：
由  并令  可得：
由  可得：
由  可得：
我们令  ，则  ，由  可得  。
将  、 、 代入  中，可得：
至此， 均表示成二阶导的关系式，将其代入  可得：
这样，我们就能构造出一个以二阶导  为未知数的线性方程组了。
自然边界：
可以看出，左侧的系数矩阵为严格对角占优矩阵，故线性方程组有唯一解。再根据  关于二阶导的表达式就能求出每个分段区间上的具体方程了。
固定边界：
由  有： ，而  ，则有：
由  有： ，将  的表达式带入可得：
则线性方程组为：
非扭结边界：
由  有： ，代入  的表达式可得：
则线性方程组为：
算法
输入： 个数据节点  
[1] 计算步长  ；
[2] 配置并求解线性方程组；
[3] 计算样条曲线的系数；
[4] 在每个分段区间  上创建三次多项式方程  。
3. 三对角线性系统的求解
对于系数矩阵为严格对角占优三对角矩阵的线性方程组，有唯一解。且其  分解、前代及回代过程，只需  次操作。
在编程方面，通常无需存储整个  系数矩阵，使用三个列向量来代替。
void tridag(int n, double a[], double b[], double c[], double r[], double u[]){
    double bet;
    double *gma = vector(1, n - 1, (double)(NULL));

    bet = b[0];
    u[0] = r[0] / bet;

    for (int i = 1; i// 分解及前代
        gma[i] = c[i - 1] / bet;
        bet = b[i] - a[i] * gma[i];
        u[i] = (r[i] - a[i] * u[i - 1]) / bet;
    }

    for (int i = n - 2; i >= 0; i--) {    // 回代
        u[i] -= gma[i + 1] * u[i + 1];
    }
}
4. 固定边界三次样条的求解
仿照前述三对角线性系统的求解方法，对于固定边界的三次样条，编写如下样条插值的功能块。
/**
 * 在处理列表函数时，程序spline仅调用一次。一旦调用过后，对任意的xx，插值函数的值通过splint获得。
 **/

/**
 * 给定数组x[0...n-1](严格单调递增)与y[0...n-1]构成一个列表函数，即y[i] = f(x[i])。
 * yp1与ypn表示插值函数分别在第1个与第n-1个点处的一阶导数。
 * 函数修改y2[0...n-1]，表示插值函数在表中点x_i处的二阶导数。
 **/
void spline(double x[], double y[], int n, double yp1, double ypn, double y2[]){
    double bet;
    double *gma = vector(1, n-1, (double)(NULL));

    bet = 2*(x[1]-x[0]);
    y2[0] = 6*((y[1]-y[0])/(x[1]-x[0]) - yp1) / bet;

    for(int i=1; i-1; i++) {    // 分解及前代
        gma[i] = (x[i]-x[i-1]) / bet;
        bet = 2*(x[i+1]-x[i-1]) - (x[i]-x[i-1]) * gma[i];
        y2[i] = (6*((y[i+1]-y[i])/(x[i+1]-x[i]) - (y[i]-y[i-1])/(x[i]-x[i-1])) - (x[i]-x[i-1]) * y2[i-1]) / bet;
    }
    gma[n-1] = (x[n-1]-x[n-2]) / bet;
    bet = 2*(x[n-1]-x[n-2]) * gma[n-1];
    y2[n-1] = (6*(ypn - (y[n-1]-y[n-2])/(x[n-1]-x[n-2])) - (x[n-1]-x[n-2]) * y2[n-2]) / bet;for(int i=n-2; i>=0; i--) {    // 回代
        y2[i] -= gma[i+1] * y2[i+1];
    }
}/**
 * 给定数组x[0...n-1]和y[0...n-1]作为列表函数，数组y2a[0...n-1]由spline程序输出给定。
 * 给定要计算的xx值，程序通过yy返回其三次样条插值结果。
 **/void splint(double x[], double y[], double y2[], int n, double xx, double *yy){int k, kl=0, kh=n-1;double a, b, c, d;// find the corresponding interval with binary searchingwhile(kh-kl > 1) {
        k = (kh+kl)/2;if (x[k] > xx) {
            kh = k;
        } else {
            kl = k;
        }
    }
    a = y[kl];
    b = (y[kh]-y[kl])/(x[kh]-x[kl]) - (x[kh]-x[kl])*y2[kl]/2 - (x[kh]-x[kl])*(y2[kh]-y2[kl])/6;
    c = y2[kl]/2;
    d = (y2[kh]-y2[kl])/(x[kh]-x[kl])/6;
    *yy = a + b*(xx-x[kl]) + c*(xx-x[kl])*(xx-x[kl]) + d*(xx-x[kl])*(xx-x[kl])*(xx-x[kl]);
}
测试
针对上述功能块  和  ，进行如下测试：在  上选取若干个控制点，计算控制点之外的插值点。
#include 
#include 
#include "nrutil.h"

int main(){
    FILE *fp;
    int n = 21;

    double xx;
    double yy;
    double yp1 = -20;
    double ypn = 20;

    double *x = vector(0, n-1, (double)(NULL));
    double *y = vector(0, n-1, (double)(NULL));
    double *y2 = vector(0, n-1, (double)(NULL));

    // select n control points from f(x)=x^2
    fp = fopen("control.dat","w");
    for (int i=0; i        x[i] = i-10.0;
        y[i] = (i-10.0)*(i-10.0);
        fprintf(fp,"%lf\t%lf\n",x[i],y[i]);
    }
    fclose(fp);

    spline(x, y, n, yp1, ypn, y2);

    // calculate the interpolating points
    fp=fopen("interpolate.dat","w");
    xx = -10.0;
    while (xx <= 10.0) {
        splint(x, y, y2, n, xx, &yy);
        fprintf(fp,"%lf\t%lf\n",xx,yy);
        xx += 0.1;
    }
    fclose(fp);

    return 0;
}
画图脚本：
set xlabel 'x'
set ylabel 'y'

set multiplot layout 1,1 title 'cubic spline'

plot "control.dat" u 1:2 t 'control' w points pointtype 7 linecolor 'red' linewidth 1 pointsize 1,\
     "interpolate.dat" u 1:2 t 'interpolatation' w points pointtype 6 linecolor 'blue' linewidth 1 pointsize 1 

unset multiplot
效果图：
5. 附录
5.1 龙格库塔现象
龙格现象是在一组等间距插值点上使用高次多项式插值时出现的区间边缘处的振荡问题，它表明使用高次多项式插值并不总是能提高准确性。随着多项式阶次的增加，以这种方式产生的多项式  实际上可能偏离  ，通常发生在靠近插值点的端点。
解决龙格现象的方法：
使用分段多项式样条可以避免龙格库塔现象。如果要减小插值误差，那么可以增加构成样条的多项式的数目，而不必是增加多项式的阶次。