Cubic Spline Interpolation

对于$n+1$给定点的数据集${\displaystyle \{x_i\}},i=0,1,\cdots ,n$,我们可以用$n$段三次多项式在数据点之间构建一个三次样条.如果$$S(x)=\{\begin{array}{c}{S}_0(x),x\in [x_0,x_1]\\ S_1(x),x\in [x_1,x_2]\\ \cdots \\ S_{n-1}(x),x\in [x_{n-1},x_n]\end{array}$$表示对函数$f$进行插值的样条函数,那么需要:

• 插值特性,${\displaystyle S(x_i)=f(x_i)},i=0,1,\cdots ,n$
• 样条相连,${\displaystyle S_{i-1}(x_i)=S_i(x_i),i=1,\dots ,n-1}$
• 二阶可导,${\displaystyle S_{i-1}^{\prime }(x_i)=S_i^{\prime }(x_i)}$以及$S_{i-1}^{″}(x_i)=S_i^{″}(x_i),i=1,\dots ,n-1$.

由于每个三次多项式需要四个条件才能确定曲线形状,所以对于组成$S$的$n$个三次多项式来说,这就意味着需要$4n$个条件才能确定这些多项式.但是,插值特性只给出了$n+1$个条件,内部数据点给出$n+1-2=n-1$个条件,总计是$4n-2$个条件.我们还需要另外两个条件,根据不同的因素我们可以使用不同的条件.
其中一项选择条件可以得到给定$u$与$v$的钳位三次样条(clamped cubic spline),$$S^{\prime }(x_0)=u$$$$S^{\prime }(x_k)=v$$另外,我们可以设$$S^{″}(x_0)=S^{″}(x_n)=0$$这样就得到自然三次样条.自然三次样条几乎等同于样条设备生成的曲线.
在这些所有的二次连续可导函数中,钳位与自然三次样条可以得到相对于待插值函数$f$的最小震荡.
如果选择另外一些条件,$$S(x_0)=S(x_n)$$$$S^{\prime }(x_0)=S^{\prime }(x_n)$$$$S^{″}(x_0)=S^{″}(x_n)$$可以得到周期性的三次样条.
如果选择,$$S(x_0)=S(x_n)$$$$S^{\prime }(x_0)=S^{\prime }(x_n)$$$$S^{″}(x_0)=f^{\prime }(x_0),S^{″}(x_n)=f^{\prime }(x_n)$$可以得到complete三次样条.

---

Algorithm to find the interpolating cubic spline
We wish to find each polynomial $q_i(x)$ given the points $(x_0,y_0)$ through $(x_n,y_n)$. To do this, we will consider just a single piece of the curve, $q(x)$, which will interpolate from $(x_1,y_1)$ to $(x_2,y_2)$. This piece will have slopes $k_1$ and $k_2$ at its endpoints. Or, more precisely,

$$q(x_1)=y_1,$$$$q(x_2)=y_2,$$$$q^{\prime }(x_1)=k_1,$$$$q^{\prime }(x_2)=k_2.$$The full equation $q(x)$ can be written in the symmetrical form$$\begin{array}{lcr}\text{(1)}& q(x)=(1-t(x))y_1+t(x)y_2+t(x)(1-t(x))((1-t(x))a+t(x)b),& \text{(1)}\end{array}$$where$$\begin{array}{lcr}\text{(2)}& t(x)=\frac{x-x_1}{x_2-x_1},& \text{(2)}\end{array}$$$$\begin{array}{lcr}\text{(3)}& a=k_1(x_2-x_1)-(y_2-y_1),& \text{(3)}\end{array}$$$$\begin{array}{lcr}\text{(4)}& b=-k_2(x_2-x_1)+(y_2-y_1).& \text{(4)}\end{array}$$

But what are $k_1$ and $k_2$? To derive these critical values, we must consider that$$q^{\prime }=\frac{dq}{dx}=\frac{dq}{dt}\frac{dt}{dx}=\frac{dq}{dt}\frac{1}{x_2-x_1}.$$It then follows that$$\begin{array}{lcr}\text{(5)}& q^{\prime }=\frac{y_2-y_1}{x_2-x_1}+(1-2t)\frac{a(1-t)+bt}{x_2-x_1}+t(1-t)\frac{b-a}{x_2-x_1},& \text{(5)}\end{array}$$$$\begin{array}{lcr}\text{(6)}& q^{″}=2\frac{b-2a+(a-b)3t}{{(x_2-x_1)}^2}.& \text{(6)}\end{array}$$Setting $t=0$ and $t=1$ respectively in equations $\text{(5)}$ and $\text{(6)}$, one gets from $\text{(2)}$ that indeed first derivatives $q^{\prime }(x_1)=k_1$ and $q^{\prime }(x_2)=k_2$, and also second derivatives$$\begin{array}{lcr}\text{(7)}& q^{″}(x_1)=2\frac{b-2a}{{(x_2-x_1)}^2},& \text{(7)}\end{array}$$$$\begin{array}{lcr}\text{(8)}& q^{″}(x_2)=2\frac{a-2b}{{(x_2-x_1)}^2}.& \text{(8)}\end{array}$$If now $(x_i,y_i),i=0,1,...,n$ are $n+1$ points, and$$\begin{array}{lcr}\text{(9)}& q_i=(1-t)y_{i-1}+t{y}_i+t(1-t)((1-t)a_i+t{b}_i),& \text{(9)}\end{array}$$where $i=1,2,...,n$, and $t=\frac{x-x_{i-1}}{x_i-x_{i-1}}$ are $n$ third-degree polynomials interpolating $y$ in the interval $x_{i-1}\le x\le x_i$ for $i=1,...,n$ such that $q_i^{\prime }(x_i)=q_{i+1}^{\prime }(x_i)$ for $i=1,...,n-1$, then the $n$ polynomials together define a differentiable function in the interval $x_0\le x\le x_n$, and$$\begin{array}{lcr}\text{(10)}& a_i=k_{i-1}(x_i-x_{i-1})-(y_i-y_{i-1}),& \text{(10)}\end{array}$$$$\begin{array}{lcr}\text{(11)}& b_i=-k_i(x_i-x_{i-1})+(y_i-y_{i-1})& \text{(11)}\end{array}$$for $i=1,...,n$, where$$\begin{array}{lcr}\text{(12)}& k_0=q_1^{\prime }(x_0),& \text{(12)}\end{array}$$$$\begin{array}{lcr}\text{(13)}& k_i=q_i^{\prime }(x_i)=q_{i+1}^{\prime }(x_i),i=1,\dots ,n-1,& \text{(13)}\end{array}$$$$\begin{array}{lcr}\text{(14)}& k_n=q_n^{\prime }(x_n).& \text{(14)}\end{array}$$If the sequence $k_0,k_1,\cdots ,k_n$ is such that, in addition, $q_i^{″}(x_i)=q_{i+1}^{″}(x_i)$ holds for $i=1,2,\cdots ,n-1$, then the resulting function will even have a continuous second derivative.
From $\text{(7)}$,$\text{(8)}$,$\text{(10)}$ and $\text{(11)}$ follows that this is the case if and only if$$\begin{array}{lcr}\text{(15)}& \frac{k_{i-1}}{x_i-x_{i-1}}+(\frac{1}{x_i-x_{i-1}}+\frac{1}{x_{i+1}-x_i})2k_i+\frac{k_{i+1}}{x_{i+1}-x_i}=3(\frac{y_i-y_{i-1}}{{(x_i-x_{i-1})}^2}+\frac{y_{i+1}-y_i}{{(x_{i+1}-x_i)}^2})& \text{(15)}\end{array}$$for $i=1,2,\cdots ,n-1$. The relations $\text{(15)}$ are $n-1$ linear equations for the $n+1$ values $k_0,k_1,...,k_n$.
For the elastic rulers being the model for the spline interpolation, one has that to the left of the left-most "knot" and to the right of the right-most "knot" the ruler can move freely and will therefore take the form of a straight line with $q^{″}=0$. As $q^{″}$ should be a continuous function of $x$, "natural splines" in addition to the $n-1$ linear equations $\text{(15)}$ should have$$q_1^{″}(x_0)=2\frac{3(y_1-y_0)-(k_1+2k_0)(x_1-x_0)}{{(x_1-x_0)}^2}=0,$$$$q_n^{″}(x_n)=-2\frac{3(y_n-y_{n-1})-(2k_n+k_{n-1})(x_n-x_{n-1})}{{(x_n-x_{n-1})}^2}=0,$$i.e. that
$$\begin{array}{lcr}\text{(16)}& \frac{2}{x_1-x_0}{k}_0+\frac{1}{x_1-x_0}{k}_1=3\frac{y_1-y_0}{(x_1-x_0{)}^2},& \text{(16)}\end{array}$$$$\begin{array}{lcr}\text{(17)}& \frac{1}{x_n-x_{n-1}}{k}_{n-1}+\frac{2}{x_n-x_{n-1}}{k}_n=3\frac{y_n-y_{n-1}}{(x_n-x_{n-1}{)}^2}.& \text{(17)}\end{array}$$Eventually, $\text{(15)}$ together with $\text{(16)}$ and $\text{(17)}$ constitute $n+1$ linear equations that uniquely define the $n+1$ parameters $k_0,k_1,\cdots ,k_n$.
There exist other end conditions, "clamped spline", which specifies the slope at the ends of the spline, and the popular "not-a-knot spline", which requires that the third derivative is also continuous at the $x_1$ and $x_{N-1}$ points.
For the "not-a-knot" spline, the additional equations will read:$$q_1^{‴}(x_1)=q_2^{‴}(x_1)\Rightarrow \frac{1}{\mathrm{\Delta }{x}_1^2}{k}_0+(\frac{1}{\mathrm{\Delta }{x}_1^2}-\frac{1}{\mathrm{\Delta }{x}_2^2})k_1-\frac{1}{\mathrm{\Delta }{x}_2^2}{k}_2=2(\frac{\mathrm{\Delta }{y}_1}{\mathrm{\Delta }{x}_1^3}-\frac{\mathrm{\Delta }{y}_2}{\mathrm{\Delta }{x}_2^3}),$$$$q_{n-1}^{‴}(x_{n-1})=q_n^{‴}(x_{n-1})\Rightarrow \frac{1}{\mathrm{\Delta }{x}_{n-1}^2}{k}_{n-2}+(\frac{1}{\mathrm{\Delta }{x}_{n-1}^2}-\frac{1}{\mathrm{\Delta }{x}_n^2})k_{n-1}-\frac{1}{\mathrm{\Delta }{x}_n^2}{k}_n=2(\frac{\mathrm{\Delta }{y}_{n-1}}{\mathrm{\Delta }{x}_{n-1}^3}-\frac{\mathrm{\Delta }{y}_n}{\mathrm{\Delta }{x}_n^3}),$$where $\mathrm{\Delta }{x}_i=x_i-x_{i-1},\mathrm{\Delta }{y}_i=y_i-y_{i-1}$.

Example
In case of three points the values for $k_0,k_1,k_2$ are found by solving the tridiagonal linear equation system$$\left[\begin{array}{ccc}{a}_{11}& a_{12}& 0\\ a_{21}& a_{22}& a_{23}\\ 0& a_{32}& a_{33}\end{array}\right]\left[\begin{array}{c}{k}_0\\ k_1\\ k_2\end{array}\right]=\left[\begin{array}{c}{b}_1\\ b_2\\ b_3\end{array}\right]$$ with
$a_{11}=\frac{2}{x_1-x_0},$ $a_{12}=\frac{1}{x_1-x_0},$ $a_{21}=\frac{1}{x_1-x_0},$ $a_{22}=2(\frac{1}{x_1-x_0}+\frac{1}{x_2-x_1}),$ $a_{23}=\frac{1}{x_2-x_1},$ $a_{32}=\frac{1}{x_2-x_1},$ $a_{33}=\frac{2}{x_2-x_1},$ $b_1=3\frac{y_1-y_0}{(x_1-x_0{)}^2},$ $b_2=3(\frac{y_1-y_0}{{(x_1-x_0)}^2}+\frac{y_2-y_1}{{(x_2-x_1)}^2}),$ $b_3=3\frac{y_2-y_1}{(x_2-x_1{)}^2}.$
For the three points
$(-1,0.5),(0,0),(3,3),$ one gets that
$k_0=-0.6875,k_1=-0.1250,k_2=1.5625,$
and from $\text{(10)}$ and $\text{(11)}$ that
$a_1=k_0(x_1-x_0)-(y_1-y_0)=-0.1875,$
$b_1=-k_1(x_1-x_0)+(y_1-y_0)=-0.3750,$
$a_2=k_1(x_2-x_1)-(y_2-y_1)=-3.3750,$
$b_2=-k_2(x_2-x_1)+(y_2-y_1)=-1.6875.$
In the figure, the spline function consisting of the two cubic polynomials $q_1(x)$ and $q_2(x)$ given by $\text{(9)}$ is displayed.