插值是离散函数逼近的重要方法，利用它可通过函数在有限个点处的取值状况，估算出函数在其他点处的近似值。与拟合不用经过每个已知点不同，插值需要经过每个已知点，另外并不是阶数越高越好，因为高阶插值容易出现龙格现象，即插值后在区间两端点处波动极大，产生明显的震荡。三次样条插值作为一种常见的插值方法，这里记录一下其基本概念及求解过程。
一、基本概念
设在区间\([a, b]\)上存在\(n+1\)个已知数据点如下，其把\([a, b]\)分成了\(n\)个子区间\([x_0, x_1], [x_1, x_2],\ ... ,\ [x_{n-1}, x_n]\)，
\[ \begin{align} x &: a \le x_0 < x_1 < ...< x_n \le b\\ f(x) = y &: \quad\quad y_0\ \quad y_1\ \quad ... \quad y_n \end{align} \]如果函数\(S(x)\)满足以下三个条件，则称\(S(x)\)为\(f(x)\)关于节点\(x_0, x_1, ..., x_n\)的三次样条函数。
- 在每个子区间上都是一个不超过3次的多项式，即 \(S_i(x_i) = a_i + b_i(x-x_i) + c_i(x-x_i)^2 + d_i(x-x_i)^3\)；
- 在整个区间\([a,b]\)上连续且光滑，即一阶导数\(S^{'}(x_i)\)和二阶导数\(S^{''}(x_i)\)存在且连续；
- 满足插值条件，即\(S_i(x_i) = y_i,\ i = 0, 1, 2,..., n\)。
二、求解过程
2.1 获取方程组
\(n\)个子区间中每一个都有四个未知数\((a_i, b_i, c_i, d_i)\)，所以共有\(4n\)个未知参数需要求解。
插值条件
曲线在整个区间\([a, b]\)上所有点都满足插值条件，所以\(S_i(x_i) = y_i,\ i = 0, 1, 2,..., n-1,\ S_{n-1}(x_n) = y_n\)，依此可得到\(n+1\)个方程。
曲线连续
曲线在整个区间\([a, b]\)上连续，表明在端点\(i = 1, 2, ... , n-1\)处两边函数值相等，即：\(S_{i-1}(x_i) = S_{i}(x_{i})\)，其等价于\(S_i(x_{i+1}) = S_{i+1}(x_{i+1}) = y_{i+1}, i = 0, 1, 2, ..., n-2\)，这里可得\(n-1\)个方程。
曲线光滑
曲线在整个区间\([a, b]\)上光滑，表明在端点\(i = 1, 2, ... , n-1\)两边的一阶导数和二阶导数存在且相等，即
\[S^{'}_{i-1}(x_i) = S^{'}_{i}(x_i), i = 1, 2, ..., n-1.\ S^{''}_{i-1}(x_i) = S^{''}_{i}(x_i), i = 1, 2, ..., n-1 \]其等价于
\[S^{'}_{i}(x_{i+1}) = S^{'}_{i+1}(x_{i+1}), i = 0, 1, ..., n-2.\ S^{''}_{i}(x_{i+1}) = S^{''}_{i+1}(x_{i+1}), i = 0, 1, ..., n-2 \]共\(2(n-1)\)个方程。
区间\([a, b]\)左右两端点特性
由2.1---2.3一共可以得到\(n+1 + n-1 + 2(n-1) = 4n - 2\)个方程，要求解4n个未知数，还需要至少两个方程，所以这里可以考虑在两端点\(i = 0, n\)处的特性，加上边界处这两个方程可以对\(4n\)个参数进行求解。一般有三种边界条件：自然边界（Natural Spline），固定边界（Clamped Spline），非节点边界（Not-A-Knot Spline）。
-
自然边界
指定端点二阶导数为0，即\(S^{''}_0(x_0) = S^{''}_{n-1}(x_n)=0\)。 -
固定边界
人为指定端点一阶导数，这里分别定为\(A\)和\(B\)，即\(S^{'}_0(x_0) = A, S^{'}_{n-1}(x_n) = B\)。 -
非节点边界
强制第一个插值点的三阶导数值等于第二个点的三阶导数值，最后第一个点的三阶导数值等于倒数第二个点的三阶导数值. 即\(S^{'''}_0(x_0) = S^{'''}_1(x_1), S^{'''}_{n-2}(x_{n-1}) = S^{'''}_{n-1}(x_{n})\)。
2.1 方程推导
由\(S_i(x) = a_i + b_i(x-x_i) + c_i(x-x_i)^2 + d_i(x-x_i)^3\)可得其各阶导数如下：
\(S^{'}_i(x) = b_i + 2c_i(x-x_i) + 3d_i(x-x_i)^2\)
\(S^{''}_i(x) = 2c_i + 6d_i(x-x_i)\)
\(S^{'''}_i(x) = 6d_i\)
1、由\(S_i(x_i) = y_i,\ i = 0, 1, 2,..., n-1\)可得：
\(S_i(x_i) = a_i + b_i(x_i-x_i) + c_i(x_i-x_i)^2 + d_i(x_i-x_i)^3 = y_i\)，化简可得：
\(a_i = y_i \ \ \ \ \ (1)\)
2、由\(S_i(x_{i+1}) = y_{i+1},\ i = 0, 1, 2, ..., n - 2\)可得：
\(S_i(x_{i+1}) = a_i + b_i(x_{i+1}-x_i) + c_i(x_{i+1}-x_i)^2 + d_i(x_{i+1}-x_i)^3 = y_{i+1}\)，令\(h_i = x_{i+1}-x_i\)，则可简写为
\(a_i + b_ih_i + c_ih^2_i + d_ih^3_i = y_{i+1} \ \ \ \ \ (2)\)
3、由\(S^{'}_{i}(x_{i+1}) = S^{'}_{i+1}(x_{i+1}), i = 0, 1, ..., n-2\)可得：
\(S^{'}_i(x_{i+1}) = b_i + 2c_i(x_{i+1}-x_i) + 3d_i(x_{i+1}-x_i)^2 = b_i + 2c_ih_i + 3d_ih^2_i\)
\(S^{'}_{i+1}(x_{i+1}) = b_{i+1} + 2c_{i+1}(x_{i+1}-x_{i+1}) + 3d_{i+1}(x_{i+1}-x_{i+1})^2 = b_{i+1}\)
所以
\(b_i + 2c_ih_i + 3d_ih^2_i - b_{i+1} = 0 \ \ \ \ \ (3)\)
4、由\(S^{''}_{i}(x_{i+1}) = S^{''}_{i+1}(x_{i+1}), i = 0, 1, ..., n-2\)可得
\(2c_i + 6d_ih_i - 2C_{i+1} = 0 \ \ \ \ \ (4)\)
通过(2)(3)(4)式可知\(b_i, c_i, d_i\)存在某种关系，为了方便求解，这里作以下变换，将\(b_i, c_i, d_i\)三个参数变换为求解\(m_i\)一个参数。
由(4)式可知：
\(d_i = \frac{2C_{i+1} - 2C_{i}}{6h_i}\), 令\(2C_i = m_i\),则：
\(d_i = \frac{m_{i+1} - m_{i}}{6h_i}\)
由(2)式可知：
\(a_i + b_ih_i + c_ih^2_i + d_ih^3_i = y_{i+1} = y_i + h_ib_i + h^2_i\frac{m_i}{2} + h^3_i\frac{m_{i+1} - m_i}{6h_i}\)，化简可得：
\(b_i = \frac{y_{i+1} - y_i}{h_i} = \frac{h_i}{2}m_i - \frac{h_i}{6}(m_{i+1} - m_i)\)
将\(b_i, c_i, d_i\)代入(3)式可知：
\(h_im_i + 2(h_i + h_{i+1})m_{i+1} + h_{i+1}m_{i+2} = 6(\frac{y_{i+2} - y_{i+1}}{h_{i+1}} - \frac{y_{i+1} - y_i}{h_i}), i = 0, 1, ..., n-2.\)
通过以上推导可知关于\(m_i\)的一个方程组，维度为\(n-1\)，再考虑边界条件便可以得到\(n+1\)个方程，进而进行求解。
-
自由边界
由\(S^{''}_0(x_0) = S^{''}_{n-1}(x_n)=0\)可知\(m_0 = m_n = 0\)，因此方程组可以写成以下形式
-
固定边界
因此方程组的系数矩阵可以改写为
-
非节点边界
因此方程组的系数矩阵可以改写为
三、算法总结
-
步骤1：分区间
将数据分成不同子区间子区间\([x_0, x_1], [x_1, x_2],\ ... ,\ [x_{n-1}, x_n]\) -
步骤2：计算步长
计算步长\(h_i = x_{i+1} - x_{i}\) -
步骤3：求解方程组获得\(m_i\)
-
步骤4：计算每个子区间的参数{a_i, b_i, c_i, d_i}
\(a_i = y_i\)
\(b_i = \frac{y_{i+1} - y_i}{h_i} = \frac{h_i}{2}m_i - \frac{h_i}{6}(m_{i+1} - m_i)\)
\(c_i = \frac{m_i}{2}\)
\(d_i = \frac{m_{i+1} - m_{i}}{6h_i}\)
四、代码实现
参考三次样条插值(Cubic Spline Interpolation)及代码实现(C语言)
参考链接
三次样条插值(Cubic Spline Interpolation)及代码实现(C语言)
三次样条（cubic spline）插值