数值计算—三次样条函数插值

一、样条函数的基本原理

已知 $n+1$ 个数据节点，分成 $n$ 个区间（Interval），在每个区间构造一个三次函数，$S_0(x)\sim S_{n-1}(x)$。

假设每一段三次函数都用如下方程来定义：
$$S_i(x)=a_i+b_i(x-x_i)+c_i(x-x_i{)}^2+d_i(x-x_i{)}^3,i=0,1,\cdots ,n-1$$

那么 $n$ 个区间对应的三次函数 $S(x)$ 的数学表达式如下：
$$S(x)=\{\begin{array}{ll}{S}_0(x)=a_0+b_0(x-x_0)+c_0(x-x_0{)}^2+d_0(x-x_0{)}^3& if{x}_0\le x\le x_1\\ S_1(x)=a_1+b_1(x-x_1)+c_1(x-x_1{)}^2+d_1(x-x_1{)}^3& if{x}_1\le x\le x_2\\ ⋮& \\ S_{n-1}(x)=a_{n-1}+b_{n-1}(x-x_{n-1})+c_{n-1}(x-x_{n-1}{)}^2+d_{n-1}(x-x_{n-1}{)}^3& if{x}_{n-1}\le x\le x_n\end{array}$$

每段三次函数 $S(x)$ 都有 $a,b,c,d$ 四个系数，上述 $n$ 段三次函数，共有 $4n$ 个系数。要想求解这 $4n$ 个系数，共需要 $4n$ 个独立的方程。下面，我们根据三次样条函数对每一段三次函数在断点处的约束，生成求解 $4n$ 个系数所需要的所有方程。

条件1：$n$ 段三次函数必须穿过所有已知节点

$$\boxed{S_i(x_i)=y_i,i=0,1,2,\cdots ,n}$$

已知 $n+1$ 个节点，总共可以生成 $n+1$ 个方程，即
$$\{\begin{array}{l}{y}_0=S_0(x_0)=a_0+b_0(x_0-x_0)+c_0(x_0-x_0{)}^2+d_0(x_0-x_0{)}^3=a_0\\ y_1=S_1(x_1)=a_1+b_1(x_1-x_1)+c_1(x_1-x_1{)}^2+d_1(x_1-x_1{)}^3=a_1\\ ⋮\\ y_{n-1}=S_{n-1}(x)=a_{n-1}+b_{n-1}(x_{n-1}-x_{n-1})+c_{n-1}(x_{n-1}-x_{n-1}{)}^2+d_{n-1}(x_{n-1}-x_{n-1}{)}^3=a_{n-1}\\ y_n=S_{n-1}(x)=a_{n-1}+b_{n-1}(x_n-x_{n-1})+c_{n-1}(x_n-x_{n-1}{)}^2+d_{n-1}(x_n-x_{n-1}{)}^3\end{array}$$

条件2：在所有节点（除了第一节点和最后一个节点）处连续

$$\boxed{S_i(x_{i+1})=S_{i+1}(x_{i+1}),i=0,1,2,\cdots ,n-2}$$

已知 $n+1$ 个节点，$n$ 个区间，因此有 $n-1$ 个衔接点，共生成 $n-1$ 个方程，即
$$\{\begin{array}{l}{S}_0(x_1)=S_1(x_1)\\ S_1(x_2)=S_2(x_2)\\ ⋮\\ S_{n-2}(x_{n-1})=S_{n-1}(x_{n-1})\end{array}$$

条件3：在所有节点（除了第一节点和最后一个节点）处可导

$$\boxed{S_i^{\prime }(x_{i+1})=S_{i+1}^{\prime }(x_{i+1}),i=0,1,2,\cdots ,n-2}$$

已知 $n+1$ 个节点，$n$ 个区间，因此有 $n-1$ 个衔接点，共生成 $n-1$ 个方程，即
$$\{\begin{array}{l}{S}_0^{\prime }(x_1)=S_1^{\prime }(x_1)\\ S_1^{\prime }(x_2)=S_2^{\prime }(x_2)\\ ⋮\\ S_{n-2}^{\prime }(x_{n-1})=S_{n-1}^{\prime }(x_{n-1})\end{array}$$

条件4：在所有节点（除了第一节点和最后一个节点）处二阶可导

$$\boxed{S_i^{\prime \prime }(x_{i+1})=S_{i+1}^{\prime \prime }(x_{i+1}),i=0,1,2,\cdots ,n-2}$$

已知 $n+1$ 个节点，$n$ 个区间，因此有 $n-1$ 个衔接点，共生成 $n-1$ 个方程，即
$$\{\begin{array}{l}{S}_0^{\prime \prime }(x_1)=S_1^{\prime \prime }(x_1)\\ S_1^{\prime \prime }(x_2)=S_2^{\prime \prime }(x_2)\\ ⋮\\ S_{n-2}^{\prime \prime }(x_{n-1})=S_{n-1}^{\prime \prime }(x_{n-1})\end{array}$$

综上，第一个约束条件得到了 $n+1$ 个方程，第二、三、四个约束条件分别得到了 $n-1$ 个方程，总计 $4n-2$ 个方程，还缺少 2 个方程。

二、求解方程组

函数 $S(x)$ 的一阶导数和二阶导数为
$$\{\begin{array}{l}{S}_i(x)=a_i+b_i(x-x_i)+c_i(x-x_i{)}^2+d_i(x-x_i{)}^3\\ S_i^{\prime }(x)=b_i+2c_i(x-x_i)+3d_i(x-x_i{)}^2\\ S_i^{\prime \prime }(x)=2c_i+6d_i(x-x_i)\end{array}$$

根据条件 1， $S_i(x_i)=y_i$，推导得出：
$$\boxed{a_i=y_i}$$

根据条件 2， $S_i(x_{i+1})=S_{i+1}(x_{i+1})$（即 $S_i(x_{i+1})=y_{i+1}$）推导得出：
$$\boxed{a_i+h_i{b}_i+h_i^2{c}_i+h_i^3{d}_i=y_{i+1}}$$

$$\begin{gathered} 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 \\ 令h_i=x_{i+1}-x_i \\ {S}_i(x_{i+1})=a_i+b_i{h}_i+c_i{h}_i^2+d_i{h}_i^3 \\ 又S_{i+1}(x_{i+1})=y_{i+1} \\ 故，a_i+b_i{h}_i+c_i{h}_i^2+d_i{h}_i^3=y_{i+1} \end{gathered}$$

根据条件 3 ，$S_i^{\prime }(x_{i+1})=S_{i+1}^{\prime }(x_{i+1})$ 推导得出：
$$\boxed{b_i+2h_i{c}_i+3h_i^2{d}_i-b_{i+1}=0}$$

$$\begin{array}{rl}{S}_i^{\prime }(x_{i+1})& =b_i+2c_i(x_{i+1}-x_i)+3d_i(x_{i+1}-x_i{)}^2\\ & =b_i+2c_i{h}_i+3d_i{h}_i^2\\ S_{i+1}^{\prime }(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_i{h}_i+3d_i{h}_i^2-b_{i+1}=0\end{array}$$

根据条件 4 ，$S_i^{\prime \prime }(x_{i+1})=S_{i+1}^{\prime \prime }(x_{i+1})$ 推导得出：
$$\boxed{2c_i+6h_i{d}_i-2c_{i+1}=0}$$

$$\begin{array}{rl}{S}_i^{\prime \prime }(x_{i+1})& =2c_i+6d_i(x_{i+1}-x_i)\\ & =2ci+6d_i{h}_i\\ S_{i+1}^{\prime \prime }(x_{i+1})& =2c_{i+1}+6d_{i+1}(x_{i+1}-x_{i+1})\\ & =2c_{i+1}\\ 故，& 2c_i+6d_i{h}_i-2c_{i+1}=0\end{array}$$

综上，有
$$\{\begin{array}{ll}{a}_i=y_i& ①\\ a_i+h_i{b}_i+h_i^2{c}_i+h_i^3{d}_i=y_{i+1}& ②\\ b_i+2h_i{c}_i+3h_i^2{d}_i-b_{i+1}=0& ③\\ 2c_i+6h_i{d}_i-2c_{i+1}=0& ④\end{array}$$

令 $m_i=S_i^{\prime \prime }(x_i)=2c_i$，可得
$$\boxed{c_i=\frac{m_i}{2}}$$

由 ④ 可得
$$\boxed{d_i=\frac{m_{i+1}-m_i}{6h_i}}$$

由 ② 可得
$$\boxed{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)}$$

由 ③ 可得
$$\boxed{h_i{m}_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}]}$$

三、端点条件（自由边界）

指定端点二阶导数为 0，即
$$S_0^{\prime \prime }(x_0)=S_{n-1}^{\prime \prime }(x_n)=0$$

可得：
$$\begin{gathered} 2c_0+6d_0(x_0-x_0)=0 \\ 2c_0=0 \\ \boxed{m_0=0} \\ 2c_{n-1}+6d_{n-1}(x_n-x_{n-1})=0 \\ 2(\frac{m_{n-1}}{2})+6(\frac{m_n-m_{n-1}}{6h_{n-1}})h_{n-1}=0 \\ {m}_{n-1}+m_n-m_{n-1}=0 \\ \boxed{m_n=0} \end{gathered}$$

综上，用矩阵形式表示为：
$$\left[\begin{array}{cccccc}1& 0& 0& 0& \cdots & 0\\ h_0& 2(h_0+h_1)& h_1& 0& \cdots & 0\\ 0& h_1& 2(h_1+h_2)& h_2& \cdots & 0\\ 0& 0& h_2& 2(h_2+h_3)& \cdots & 0\\ ⋮& ⋮& ⋮& ⋮& & ⋮\\ 0& 0& 0& 0& \cdots & 1\end{array}\right]\left[\begin{array}{c}{m}_0\\ m_1\\ m_2\\ m_3\\ ⋮\\ m_n\end{array}\right]=6\left[\begin{array}{c}0\\ \frac{y_2-y_1}{h_1}-\frac{y_1-y_0}{h_0}\\ \frac{y_3-y_2}{h_2}-\frac{y_2-y_1}{h_1}\\ \frac{y_4-y_3}{h_3}-\frac{y_3-y_2}{h_2}\\ ⋮\\ 0\end{array}\right]$$

python 代码实现

import numpy as np
import matplotlib.pyplot as plt


def cubic_spline_interpolation(x, y):
    n = len(x)  # n 个插值节点，将曲线分成 n-1 个区间，每个区间构造一个三次函数（4个参数），则共有 4（n-1）个参数待求解
    h = np.zeros(n - 1)
    for i in range(0, n - 1):
        h[i] = x[i + 1] - x[i]
    A = np.zeros((n, n))
    B = np.zeros(n)
    for i in range(1, n - 1):
        A[i, i - 1] = h[i - 1]
        A[i, i] = 2 * (h[i - 1] + h[i])
        A[i, i + 1] = h[i]
        B[i] = 6 * ((y[i + 1] - y[i]) / h[i] - (y[i] - y[i - 1]) / h[i - 1])
    A[0, 0] = A[n - 1, n - 1] = 1
    B[0] = B[n - 1] = 0
    m = np.linalg.inv(A) @ B.T

    b = np.zeros(n - 1)
    d = np.zeros(n - 1)
    a = y
    c = m / 2
    for i in range(0, n - 1):
        d[i] = (m[i + 1] - m[i]) / (6 * h[i])
        b[i] = (y[i + 1] - y[i]) / h[i] - h[i] * m[i] / 2 - h[i] * (m[i + 1] - m[i]) / 6

    return a, b, c, d


x = np.array([1, 2, 3, 4, 5])
y = np.array([6, 3, 7, 2, 5])

a, b, c, d = cubic_spline_interpolation(x, y)

x_new = np.linspace(1, 5, 100)

num = len(x_new)
y_new = np.zeros(num)

for i, xi in enumerate(x_new):
    if x[-1] >= xi >= x[0]:
        j = np.searchsorted(x, xi)  # 得到 xi 在数组 x 中哪个区间中
        if j != 0:
            dx = x[j] - x[j - 1]
            dy = y[j] - y[j - 1]
            y_new[i] = a[j - 1] + b[j - 1] * (xi - x[j - 1]) + c[j - 1] * (xi - x[j - 1]) ** 2 + d[j - 1] * (
                    xi - x[j - 1]) ** 3
        else:
            y_new[i] = a[0]
    elif xi < x[0]:
        y_new[i] = a[0] + b[0] * (xi - x[0]) + c[0] * (xi - x[0]) ** 2 + d[0] * (xi - x[0]) ** 3
    elif xi > x[-1]:
        y_new[i] = a[-1] + b[-1] * (xi - x[-1]) + c[-1] * (xi - x[-1]) ** 2 + d[-1] * (xi - x[-1]) ** 3

plt.plot(x, y, 'ro', label='Original Data')
plt.plot(x_new, y_new, label='Spline Interpolation')
plt.legend()
plt.show()

结果显示