原文参考 Cubic Spline Interpolation

样条来源于早期工程制图，为了将一些固定点连成一条光滑的曲线，采用具有弹性的木条固定在这些点上，如图1所示，通过样条画出来的曲线不仅经过各固定点，而且连续光滑。后来这一技术发展成了数学工具——样条函数，在很多领域都有应用，样条插值就是典型的应用场景。 图1，利用样条进行工程绘图

介绍原理前先了解一下样条函数，顾名思义就是描述图1中那个木条的函数，在数学上把它定义成一个分段多项式函数，以三次多项式为例：

f ( x ) = { a 1 x 3 + b 1 x 2 + c 1 x + d 1 if  x ∈ [ x 1 , x 2 ] a 2 x 3 + b 2 x 2 + c 2 x + d 2 if  x ∈ ( x 2 , x 3 ] … a n x 3 + b n x 2 + c n x + d n if  x ∈ ( x n , x n + 1 ]   . (1) f\left(x\right)=\begin{cases} a_1x^3+b_1x^2+c_1x+d_1&\text{if }x\in\left[x_1,x_2\right]\\ a_2x^3+b_2x^2+c_2x+d_2&\text{if }x\in\left(x_2,x_3\right]\\ \dots\\ a_nx^3+b_nx^2+c_nx+d_n&\text{if }x\in\left(x_{n},x_{n+1}\right]\,. \end{cases}\tag{1} f(x)=⎩⎪⎪⎪⎪⎨⎪⎪⎪⎪⎧​a1​x3+b1​x2+c1​x+d1​a2​x3+b2​x2+c2​x+d2​…an​x3+bn​x2+cn​x+dn​​if x∈[x1​,x2​]if x∈(x2​,x3​]if x∈(xn​,xn+1​].​(1)

也就是说，每两个点之间用一个多项式来表示，这些多项式的阶数相同，但是系数是不一样的。样条插值的目的就是求这些多项式系数，已知系数，那些中间点的值就能唯一确定了。不难想象，随着多项式阶数的提升，端点处会更加平滑，但计算量也越大。 图2，不同阶数样条插值结果

简单来说样条插值的原理遵循四个字：连续、光滑。还有一个关键要素：边界条件。以三次样条插值（Cubic Spline Interpolation）为例，它的样条函数如公式(1)所示， ( x 1 , y 1 ) , . . . ( x n + 1 , y n + 1 ) (x_1,y_1),...(x_{n+1},y_{n+1}) (x1​,y1​),...(xn+1​,yn+1​)是已知的n+1个数据点，整个函数有 4 n 4n 4n个未知数，既 a 1 , b 1 , c 1 , d 1 , . . . , a n , b n , c n , d n a_1,b_1,c_1,d_1,...,a_n,b_n,c_n,d_n a1​,b1​,c1​,d1​,...,an​,bn​,cn​,dn​。若想求得 4 n 4n 4n个未知系数，需要构建 4 n 4n 4n个方程。 首先，根据连续性原则，每个分段函数都经过其两侧端点，因此，可得到 2 n 2n 2n个方程：

a 1 x 1 3 + b 1 x 1 2 + c 1 x 1 + d 1 = y 1 a 1 x 2 3 + b 1 x 2 2 + c 1 x 2 + d 1 = y 2 a 2 x 2 3 + b 2 x 2 2 + c 2 x 2 + d 2 = y 2 a 2 x 3 3 + b 2 x 3 2 + c 2 x 3 + d 2 = y 3 … a n x n 3 + b n x n 2 + c n x n + d n = y n a n x n + 1 3 + b n x n + 1 2 + c n x n + 1 + d n = y n + 1   . (2) a_1x_1^3+b_1x_1^2+c_1x_1+d_1=y_1\\ a_1x_2^3+b_1x_2^2+c_1x_2+d_1=y_2\\ a_2x_2^3+b_2x_2^2+c_2x_2+d_2=y_2\\ a_2x_3^3+b_2x_3^2+c_2x_3+d_2=y_3\\ \dots\\ a_nx_n^3+b_nx_n^2+c_nx_n+d_n=y_n\\ a_nx_{n+1}^3+b_nx_{n+1}^2+c_nx_{n+1}+d_n=y_{n+1}\,. \tag{2} a1​x13​+b1​x12​+c1​x1​+d1​=y1​a1​x23​+b1​x22​+c1​x2​+d1​=y2​a2​x23​+b2​x22​+c2​x2​+d2​=y2​a2​x33​+b2​x32​+c2​x3​+d2​=y3​…an​xn3​+bn​xn2​+cn​xn​+dn​=yn​an​xn+13​+bn​xn+12​+cn​xn+1​+dn​=yn+1​.(2)

其次，根据光滑行原则，相邻的两个分段函数连接处低阶导数相等，在三次样条函数中，它的低阶导数有一阶导数： d d x f i ( x ) = 3 a i x 2 + 2 b i x + c i \frac{d}{dx}f_i\left(x\right)=3a_ix^2+2b_ix+c_i dxd​fi​(x)=3ai​x2+2bi​x+ci​和二阶导数 d 2 d x 2 f i ( x ) = 6 a i x + 2 b i \frac{d^2}{dx^2}f_i\left(x\right)=6a_ix+2b_i dx2d2​fi​(x)=6ai​x+2bi​，这样，又得到 2 ( n − 1 ) 2(n-1) 2(n−1)个方程：

3 a 1 x 2 2 + 2 b 1 x 2 + c 1 = 3 a 2 x 2 2 + 2 b 2 x 2 + c 2 3 a 2 x 3 2 + 2 b 2 x 3 + c 2 = 3 a 3 x 3 2 + 2 b 3 x 3 + c 3 … 3 a n − 1 x n 2 + 2 b n − 1 x n + c n − 1 = 3 a n x n 2 + 2 b n x n + c n   . (3) 3a_1x_2^2+2b_1x_2+c_1=3a_2x_2^2+2b_2x_2+c_2\\ 3a_2x_3^2+2b_2x_3+c_2=3a_3x_3^2+2b_3x_3+c_3\\ \dots\\ 3a_{n-1}x_n^2+2b_{n-1}x_n+c_{n-1}=3a_nx_n^2+2b_nx_n+c_n\,.\tag{3} 3a1​x22​+2b1​x2​+c1​=3a2​x22​+2b2​x2​+c2​3a2​x32​+2b2​x3​+c2​=3a3​x32​+2b3​x3​+c3​…3an−1​xn2​+2bn−1​xn​+cn−1​=3an​xn2​+2bn​xn​+cn​.(3)

6 a 1 x 2 + 2 b 1 = 6 a 2 x 2 + 2 b 2 6 a 2 x 3 + 2 b 2 = 6 a 3 x 3 + 2 b 3 … 6 a n − 1 x n + 2 b n − 1 = 6 a n x n + 2 b n   . (4) 6a_1x_2+2b_1=6a_2x_2+2b_2\\ 6a_2x_3+2b_2=6a_3x_3+2b_3\\ \dots\\ 6a_{n-1}x_n+2b_{n-1}=6a_nx_n+2b_n\,.\tag{4} 6a1​x2​+2b1​=6a2​x2​+2b2​6a2​x3​+2b2​=6a3​x3​+2b3​…6an−1​xn​+2bn−1​=6an​xn​+2bn​.(4)

通过上面的连续、光滑原则得到的公式(2)、(3)、(4)总共包含 4 n − 2 4n-2 4n−2个方程，还不足以得到 4 n 4n 4n个系数的解。由于在公式(3)(4)中，第1个点和第n+1个点处在边界，只被使用了一次，导致整个方程组还差2个方程才能求解。此时，可以假设不同的边界条件来构成 4 n 4n 4n个方程。需要注意的是，不同的边界条件导致插值结果也有所不同。

Natural Spline 假设第1个和最后一个端点的二阶导为0，即：

6 a 1 x 1 + 2 b 1 = 0 6 a n x n + 1 + 2 b n = 0   . 6a_1x_1+2b_1=0\\ 6a_nx_{n+1}+2b_n=0\,.\\ 6a1​x1​+2b1​=06an​xn+1​+2bn​=0.

Not-a-Knot Spline 假设第一段和第二段函数的三阶导在端点处相等；第n-1段和第n段函数的三阶导在端点处也相等，既：

d 3 d x 3 f 1 ( x ) = d 3 d x 3 f 2 ( x ) ∣ x = x 2 d 3 d x 3 f n − 1 ( x ) = d 3 d x 3 f n ( x ) ∣ x = x n   . \frac{d^3}{dx^3}f_1\left(x\right)=\frac{d^3}{dx^3}f_2\left(x\right)\vert_{x=x_2}\\ \frac{d^3}{dx^3}f_{n-1}\left(x\right)=\frac{d^3}{dx^3}f_n\left(x\right)\vert_{x=x_n}\,. dx3d3​f1​(x)=dx3d3​f2​(x)∣x=x2​​dx3d3​fn−1​(x)=dx3d3​fn​(x)∣x=xn​​.

推到后等价于：

a 1 = a 2 a n − 1 = a n   ; a_1=a_2\\ a_{n-1}=a_n\,; a1​=a2​an−1​=an​;

Periodic Spline 假设最后一个分段函数的一阶段和二阶导于第一个分段的一阶导和二阶导相等，这种假设特别适合用在周期函数中，尤其是第一点与最后一点相等时

3 a 1 x 1 2 + 2 b 1 x 1 + c 1 = 3 a n x n + 1 2 + 2 b n x n + 1 + c n 6 a 1 x 1 + 2 b 1 = 6 a n x n + 1 + 2 b n   . 3a_1x_1^2+2b_1x_1+c_1=3a_nx_{n+1}^2+2b_nx_{n+1}+c_n\\ 6a_1x_1+2b_1=6a_nx_{n+1}+2b_n\,.\\ 3a1​x12​+2b1​x1​+c1​=3an​xn+12​+2bn​xn+1​+cn​6a1​x1​+2b1​=6an​xn+1​+2bn​.

Quadratic Spline 是最简单的假设，强制第一个与最后一个分段函数的第一个系数为0，即：

a 1 = a n = 0 a_1=a_n=0 a1​=an​=0

给定5个数据点 D = { ( 0 , 21 ) , ( 1 , 24 ) , ( 2 , 24 ) , ( 3 , 18 ) , ( 4 , 16 ) }   , \mathcal{D}=\left\{\left(0,21\right),\left(1,24\right),\left(2,24\right),\left(3,18\right),\left(4,16\right)\right\}\,, D={(0,21),(1,24),(2,24),(3,18),(4,16)},利用自然边界(Natural Spline)条件，求三次样条插值函数。5个数据点包含4个区间，所以样条函数是4分段函数，设每个分段函数表达式为： f i = a i x 3 + b i x 2 + c i x + d i , i = 1 , 2 , 3 , 4 f_i=a_ix^3+b_ix^2+c_ix+d_i,i=1,2,3,4 fi​=ai​x3+bi​x2+ci​x+di​,i=1,2,3,4，根据连续规则，得到8个方程：

f 1 ( 0 ) = 21 f 1 ( 1 ) = 24 f 2 ( 1 ) = 24 f 2 ( 2 ) = 24 f 3 ( 2 ) = 24 f 3 ( 3 ) = 18 f 4 ( 3 ) = 18 f 4 ( 4 ) = 16   . f_1\left(0\right)=21\\ f_1\left(1\right)=24\\ f_2\left(1\right)=24\\ f_2\left(2\right)=24\\ f_3\left(2\right)=24\\ f_3\left(3\right)=18\\ f_4\left(3\right)=18\\ f_4\left(4\right)=16\,. f1​(0)=21f1​(1)=24f2​(1)=24f2​(2)=24f3​(2)=24f3​(3)=18f4​(3)=18f4​(4)=16.

根据光滑原则，各分段函数一阶导与二阶导在端点处相等得到6个方程：

3 a 1 + 2 b 1 + c 1 = 3 a 2 + 2 b 2 + c 2 12 a 2 + 8 b 2 + c 2 = 12 a 3 + 8 b 3 + c 3 27 a 3 + 18 b 3 + c 3 = 27 a 4 + 18 b 4 + c 4 6 a 1 + 2 = 6 a 2 + 2 12 a 2 + 2 = 12 a 3 + 2 18 a 3 + 2 = 18 a 4 + 2 3a_1+2b_1+c_1=3a_2+2b_2+c_2\\ 12a_2+8b_2+c_2=12a_3+8b_3+c_3\\ 27a_3+18b_3+c_3=27a_4+18b_4+c_4\\ 6a_1+2=6a_2+2\\ 12a_2+2=12a_3+2\\ 18a_3+2=18a_4+2\\ 3a1​+2b1​+c1​=3a2​+2b2​+c2​12a2​+8b2​+c2​=12a3​+8b3​+c3​27a3​+18b3​+c3​=27a4​+18b4​+c4​6a1​+2=6a2​+212a2​+2=12a3​+218a3​+2=18a4​+2

根据Natural Spline假设条件，得到2个方程：

2 b 1 = 0 24 a 4 + 2 b 4 = 0   . 2b_1=0\\ 24a_4+2b_4=0\,. 2b1​=024a4​+2b4​=0.

以上方程可以联立成为一个包含16个方程得方程组，该方程组可进一步写成矩阵形式：

[ 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 8 4 2 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 8 4 2 1 0 0 0 0 0 0 0 0 0 0 0 0 27 9 3 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 27 9 3 1 0 0 0 0 0 0 0 0 0 0 0 0 64 16 4 1 3 2 1 0 − 3 − 2 − 1 0 0 0 0 0 0 0 0 0 0 0 0 0 12 4 1 0 − 12 − 4 − 1 0 0 0 0 0 0 0 0 0 0 0 0 0 27 6 1 0 − 27 − 6 − 1 0 6 2 0 0 − 6 − 2 0 0 0 0 0 0 0 0 0 0 0 0 0 0 12 2 0 0 − 12 − 2 0 0 0 0 0 0 0 0 0 0 0 0 0 0 18 2 0 0 − 18 − 2 0 0 0 2 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 24 2 0 0 ] [ a 1 b 1 c 1 d 1 a 2 b 2 c 2 d 2 a 3 b 3 c 3 d 3 a 4 b 4 c 4 d 4 ] = [ 21 24 24 24 24 18 18 16 0 0 0 0 0 0 0 0 ] \begin{bmatrix} 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 1 & 1 & 1 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 1 & 1 & 1 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 8 & 4 & 2 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 8 & 4 & 2 & 1 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 27 & 9 & 3 & 1 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 27 & 9 & 3 & 1 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 64 & 16 & 4 & 1 \\\hdashline 3 & 2 & 1 & 0 & -3 & -2 & -1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 12 & 4 & 1 & 0 & -12 & -4 & -1 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 27 & 6 & 1 & 0 & -27 & -6 & -1 & 0 \\\hdashline 6 & 2 & 0 & 0 & -6 & -2 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 12 & 2 & 0 & 0 & -12 & -2 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 18 & 2 & 0 & 0 & -18 & -2 & 0 & 0 \\\hdashline 0 & 2 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 24 & 2 & 0 & 0 \\ \end{bmatrix} \begin{bmatrix} a_1\\ b_1\\ c_1\\ d_1\\ a_2\\ b_2\\ c_2\\ d_2\\ a_3\\ b_3\\ c_3\\ d_3\\ a_4\\ b_4\\ c_4\\ d_4\\ \end{bmatrix} = \begin{bmatrix} 21\\24\\24\\24\\24\\18\\18\\16\\0\\0\\0\\0\\0\\0\\0\\0 \end{bmatrix} ⎣⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎡​0100000030060000​0100000020020020​0100000010000000​1100000000000000​00180000−3120−612000​00140000−240−22000​00120000−11000000​0011000000000000​0000827000−12270−121800​000049000−460−2200​000023000−1100000​0000110000000000​000000276400−2700−18024​00000091600−600−202​0000003400−100000​0000001100000000​​⎦⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎤​⎣⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎡​a1​b1​c1​d1​a2​b2​c2​d2​a3​b3​c3​d3​a4​b4​c4​d4​​⎦⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎤​=⎣⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎡​212424242418181600000000​⎦⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎤​

该线性系统可以利用常规方法进行求解，可得到样条函数结果：

f ( x ) = { − 0.30357 x 3 + 3.3036 x + 21 if  x ∈ [ 0 , 1 ] − 1.4821 x 3 + 3.5357 x 2 − 0.23214 x + 22.179 if  x ∈ ( 1 , 2 ] 3.2321 x 3 − 24.750 x 2 + 56.339 x − 15.536 if  x ∈ ( 2 , 3 ] − 1.4464 x 3 + 17.357 x 2 − 69.982 x + 110.79 if  x ∈ ( 3 , 4 ] . f(x) = \begin{cases}-0.30357 x^3 + 3.3036x + 21& \text{if } x \in [0,1]\\-1.4821x^3 + 3.5357x^2 -0.23214 x + 22.179& \text{if } x \in (1,2]\\3.2321x^3 -24.750x^2 + 56.339 x -15.536& \text{if } x \in (2,3]\\-1.4464x^3 + 17.357 x^2 -69.982 x + 110.79& \text{if } x \in (3,4].\end{cases} f(x)=⎩⎪⎪⎪⎪⎨⎪⎪⎪⎪⎧​−0.30357x3+3.3036x+21−1.4821x3+3.5357x2−0.23214x+22.1793.2321x3−24.750x2+56.339x−15.536−1.4464x3+17.357x2−69.982x+110.79​if x∈[0,1]if x∈(1,2]if x∈(2,3]if x∈(3,4].​

插值结果如下图所示

样条插值可以看作是一个分段函数估计问题，主要的原理就是连续性和光滑性假设，连续性在数学上的表现就是分段函数相接点处相等，光滑性表现为相接点的低阶导数相等；但光有这些条件还不够组成正定方程组，还要引入边界条件，不同的边界条件导致插值结果也有不同；此外，样条函数的阶数也是权衡计算效率和平滑性的重要参数，实际使用时根据需要来选择。