三次样条插值

上一节介绍的分段三次Hermite多项式插值在数学原理上很直接，但很多实际情况下没法应用， 因为我们可能很难甚至无法直接获得采样点的一阶导数。在 Fritsch & Carlson 发明保形特性的 PCHIP 算法之前，工业界使用更多的是一种三次样条插值(Cubic Spline Interpolation)的方法。

三次样条插值是一种典型的由工业驱动的数学建模。样条(Spline) 这个词最初来自工程制图，用于描述造船、飞机外形曲线。 为了画出平滑的曲线，工程师会将一些细木条，用砝码固定在某些特定的点上。木条在自身弹性作用下，形成的自然弯曲曲线，就是样条的物理原型。 数学上，三次样条插值就是这种物理现象的完美模拟。

三次样条函数是 Schoenberg 在 1946 年系统引入的概念。他利用了材料力学中的弹性梁理论，发现固定弹性木条在特定点上自然弯曲时，它内部所储存的弯曲能量是最小的。 数学上，所有穿过给定点的二阶连续可微的函数中，三次样条的二阶导数的平方积分最小。这在物理与数学之间建立了美妙的联系，再加上极强的应用需求，以至于三次样条插值在工业中无处不在。 典型的应用就是数控机床、机器人、自动驾驶的轨迹规划，三次样条曲线在速度和加速度上都是连续的，可以提高加工质量、机器人和车辆运动的安全性和舒适性。

1. 三次样条的数学描述

我们用 \(x_0 < x_1 < \cdots < x_{n-1} < x_n\) 共 \(n+1\) 个点上对函数 \(f(x)\) 采样。这 \(n+1\) 个点将闭区间 \([x_0, x_n]\) 分成 \(n\) 个小区间。 对于其中某个特定的子区间 \([x_i, x_{i+1}]\)，都有一个三次多项式: $$ \begin{equation}\label{f1} S_i(x) = a_i + b_i(x - x_i) + c_i(x - x_i)^2 + d_i(x - x_i)^3 \end{equation} $$ 我们需要求解出每个子区间的 4 个未知系数 \(a_i, b_i, c_i, d_i\)。需要建立如下约束:

1. 区间两端都是采样点: 对于所有 \(i = 0, 1, \cdots, n-1\)，都有 \(S_i(x_i) = f(x_i)\), \(S_i(x_{i+1}) = f(x_{i+1})\)
2. 函数值连续: 对于所有 \(i = 0, 1, \cdots, n-2\)，都有 \(S_i(x_{i+1}) = S_{i+1}(x_{i+1})\)。这点可以由 (a) 推出。
3. 一阶导数连续: 对于所有 \(i = 0, 1, \cdots, n-2\)，都有 \(S_i'(x_{i+1}) = S_{i+1}'(x_{i+1})\)。
4. 二阶导数连续: 对于所有 \(i = 0, 1, \cdots, n-2\)，都有 \(S_i''(x_{i+1}) = S_{i+1}''(x_{i+1})\)。

上述一共有 \(4n - 2\) 个约束，要求解 \(4n\) 个未知系数，还需要两个边界条件。常用的有自然样条(Natural Spline)和固定样条(Clamped Spline)两种:

1. 自然样条(Natural Spline): 令两端采样点的二阶导数为 0，即 \(S_0''(x_0) = 0, S_{n-1}''(x_n) = 0\)。其物理含义就好像，木条的两端没有任何外力作用。
2. 固定样条(Clamped Spline): 直接指定两端的一阶导数，即 \(S_0'(x_0) = f'(x_0), S_{n-1}'(x_n) = f'(x_n)\)。其物理含义是，木条两端被外力限制以固定的角度入场和出场。

从上述表述中，可以看出三次样条是 \(C^2\) 连续的，并且不需要提供采样点的一阶导数。整个 \(n\) 段三次多项式的系数，都可以通过求解线性方程组得到。 本质上它仍然是一种分段的三次 Hermite 插值多项式，只是一阶导数是由这些约束求解出来的。

2. 构造三次样条

给定采样点 \(x_0 < x_1 < \cdots < x_{n-1} < x_n\) 及其函数值 \(f(x_i), i \in \{0,1,\cdots, n\}\)， 我们可以写出 \(n\) 个式\((\ref{f1})\)中的三次多项式 \(S_i(x), i \in \{0,1,\cdots,n-1\}\)，一共有 \(4n\) 个未知系数。 根据上述 \(4n\) 个约束条件，我们可以将它们一一求解出来。

首先根据条件(a)，我们可以直接写出 \(a_i = f(x_i)\)，因为:

$$ 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 = a_i $$

条件(b) 要求采样点处函数值都是连续的，我们可以写出如下方程

$$ a_{i+1} = S_{i+1}(x_{i+1}) = 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\)，上式可以写为

$$ \begin{equation}\label{f2} a_{i+1} = a_i + b_i h_i + c_i h_i^2 + d_i h_i^3 \end{equation} $$

由于 \(S_i(x)\) 的一阶导函数为 \(S_i'(x) = b_i + 2c_i(x - x_i) + 3d_i(x - x_i)^2\)，可以直接得出 \(b_i = S_i'(x_i)\)。根据条件(c) 有:

$$ \begin{equation}\label{f3} b_{i+1} = b_i + 2c_i h_i + 3d_i h_i^2 \end{equation} $$

因为 \(S_i(x)\) 的二阶导函数为 \(S_i''(x) = 2c_i + 6d_i(x - x_i)\)，有 \(c_i = 0.5 S_i''(x_i)\)。根据条件(d) 有:

$$ \begin{equation}\label{f4} c_{i+1} = c_i + 3d_i h_i \Longrightarrow d_i = (c_{i+1} - c_i) / 3h_i \end{equation} $$

将上式中的 \(d_i\) 代入式(\(\ref{f2}\)) 有:

$$ \begin{equation}\label{f5} a_{i+1} = a_i + b_i h_i + \frac{1}{3} (c_{i+1} + 2c_i) h_i^2 \Longrightarrow \begin{cases} b_i & = \frac{1}{h_i}(a_{i+1} - a_i) - \frac{h_i}{3}(c_{i+1} + 2c_i) \\ b_{i-1} & = \frac{1}{h_{i-1}}(a_i - a_{i-1}) - \frac{h_{i-1}}{3}(c_i + 2c_{i-1}) \end{cases} \end{equation} $$

将 \(d_i\) 代入式(\(\ref{f3}\)) 有:

$$ \begin{equation}\label{f6} b_{i+1} = b_i + (c_i + c_{i+1})h_i \end{equation} $$

结合式(\(\ref{f5}\))中 \(b_i, b_{i-1}\) 的关系有 \(b_i = b_{i-1} + (c_{i-1} + c_i)h_{i-1}\)，展开有:

$$ c_{i-1}h_{i-1} + 2(h_i+h_{i-1})c_i + h_i c_{i+1} = \frac{3}{h_i}(a_{i+1} - a_i) - \frac{3}{h_{i-1}}(a_i - a_{i-1}) $$

写成向量的形式有:

$$ \begin{equation}\label{f7} \begin{bmatrix} h_{i-1} & 2(h_{i-1} + h_i) & h_i \end{bmatrix} \begin{bmatrix} c_{i-1} \\ c_i \\ c_{i+1} \end{bmatrix} = \frac{3}{h_i}(a_{i+1} - a_i) - \frac{3}{h_{i-1}}(a_i - a_{i-1}) \end{equation} $$

上式中只有 \(c_{i}, i \in \{0,\cdots,n-1\}\) 是未知的，我们一共可以写出 \(n-2\) 个上述方程，整理成矩阵的形式有:

$$ \begin{equation}\label{f8} \underbrace{\begin{bmatrix} h_0 & 2(h_0 + h_1) & h_1 & 0 & \cdots \\ 0 & h_1 & 2(h_1 + h_2) & h_2 & \\ \vdots & & \ddots & \ddots & \\ 0 & & h_{n-3} & 2(h_{n-3} + h_{n-2}) & h_{n-2} \\ \end{bmatrix}}_{\boldsymbol{H}_{(n-2)\times(n)}} \underbrace{\begin{bmatrix} c_0 \\ c_1 \\ \vdots \\ c_{n-1} \end{bmatrix}}_{\boldsymbol{C}_{n \times 1}} = \underbrace{\begin{bmatrix} 3(a_2 - a_1)/h_1 - 3(a_1 - a_0)/h_0 \\ 3(a_3 - a_2)/h_2 - 3(a_2 - a_1)/h_1 \\ \vdots \\ 3(a_{n-1} - a_{n-2})/h_{n-2} - 3(a_{n-2} - a_{n-3})/h_{n-3} \\ \end{bmatrix}}_{\boldsymbol{\beta}_{(n-2) \times 1}} \end{equation} $$

如果我们能够解出未知数\(c_{i}, i \in \{0,\cdots,n-1\}\)并代入式(\(\ref{f4},\ref{f6}\))，就可以解出所有三次多项式系数了。 显然 \(n-2\) 个方程不足以解出 \(n\) 个未知数，我们还需要附带上边界条件的约束。

3. 自然样条边界条件

我们假设两端采样点 \(x_0, x_n\) 的二阶导数为 0，即 \(f''(x_0) = 0, f''(x_n) = 0\)。那么根据 \(S_i(x)\) 的二阶导函数 \(S_i''(x) = 2c_i + 6d_i(x - x_i)\)， 我们可以轻松写出 \(c_0 = S_0''(x_0) = 0\)。此时如果曲线在 \(x >= x_n\) 处自然延伸，我们还可以写出一个三次多项式:

$$ S_n(x) = a_n + b_n(x - x_n) + c_n(x - x_n)^2 + d_n(x - x_n)^3 $$

\(S_n(x)\) 在 \(x_n\) 处保持 \(C^2\) 连续，有:

$$ \begin{cases} S_n(x_n) & = a_n = f(x_n) \\ S_n'(x_n) & = b_n \\ S_n''(x_n) & = 2c_n = 0 \\ \end{cases} $$

因此，我们可以将式(\(\ref{f8}\))扩展为:

$$ \begin{equation}\label{f9} \underbrace{\begin{bmatrix} 1 & & & & & \\ h_0 & 2(h_0 + h_1) & h_1 & 0 & \cdots & \\ 0 & h_1 & 2(h_1 + h_2) & h_2 & & \\ \vdots & & \ddots & \ddots & \ddots & \\ & & h_{n-3} & 2(h_{n-3} + h_{n-2}) & h_{n-2} & \\ & & & h_{n-2} & 2(h_{n-2} + h_{n-1}) & h_{n-1} \\ & & & & & 1 \\ \end{bmatrix}}_{\boldsymbol{H}_{(n+1)\times(n+1)}} \underbrace{\begin{bmatrix} c_0 \\ c_1 \\ c_2 \\ \vdots \\ c_{n-2} \\ c_{n-1} \\ c_n \end{bmatrix}}_{\boldsymbol{C}_{(n+1) \times 1}} = \underbrace{\begin{bmatrix} 0 \\ 3(a_2 - a_1)/h_1 - 3(a_1 - a_0)/h_0 \\ 3(a_3 - a_2)/h_2 - 3(a_2 - a_1)/h_1 \\ \vdots \\ 3(a_{n-1} - a_{n-2})/h_{n-2} - 3(a_{n-2} - a_{n-3})/h_{n-3} \\ 3(a_n - a_{n-1})/h_{n-1} - 3(a_{n-1} - a_{n-2})/h_{n-2} \\ 0 \end{bmatrix}}_{\boldsymbol{\beta}_{(n+1) \times 1}} \end{equation} $$

这样 \(\boldsymbol{H}\) 就是一个 \((n+1) \times (n+1)\) 的方阵，我们可以通过高斯消元或者 LU 分解求解线性方程组， 解出 \(c_i, i \in \{0, \cdots, n-1\}\)。我们在类 CubicSpline 中通过成员函数 InitNatureSpline 完成自然三次样条的构造。

4. 固定样条边界条件

假设采样点 \(x_0, x_n\) 的一阶导数分别为 \(f'(x_0), f'(x_n)\)。根据 \(S_i(x)\) 的一阶导函数 \(S_i'(x) = b_i + 2c_i(x - x_i) + 3d_i(x - x_i)^2\)， 有 \(b_0 = S_0'(x_0) = f'(x_0)\)。结合式(\(\ref{f5}\))，有:

$$ \begin{equation}\label{f10} b_0 = \frac{1}{h_0}(a_1 - a_0) - \frac{h_0}{3}(c_1 + 2c_0) = f'(x_0) \Longrightarrow 2h_0c_0 + h_0c_1 = \frac{3}{h_0}(a_1 - a_0) - 3 f'(x_0) \end{equation} $$

类似的，曲线在 \(x >= x_n\) 处自然延伸，有多项式 \(S_n(x) = a_n + b_n(x - x_n) + c_n(x - x_n)^2 + d_n(x - x_n)^3\)，它在 \(x_n\) 处保持 \(C^2\) 连续，有:

$$ \begin{cases} S_n(x_n) & = a_n = f(x_n) \\ S_n'(x_n) & = b_n = f'(x_n) \\ S_n''(x_n) & = 2c_n \\ \end{cases} $$

结合式(\(\ref{f5},\ref{f6}\))，有:

$$ b_n = b_{n-1} + (c_{n-1} + c_n)h_{n-1} = f'(x_n) = \frac{1}{h_{n-1}}(a_n - a_{n-1}) - \frac{h_{n-1}}{3}(c_n + 2c_{n-1}) + (c_{n-1} + c_n)h_{n-1} $$

将上式展开，有:

$$ \begin{equation}\label{f11} h_{n-1} c_{n-1} + 2h_{n-1}c_n = 3f'(x_n) - \frac{3}{h_{n-1}}(a_n - a_{n-1}) \end{equation} $$

因此，我们可以将式(\(\ref{f8}\))扩展为:

$$ \begin{equation}\label{f12} \underbrace{\begin{bmatrix} 2 h_0 & h_0 & & & & \\ h_0 & 2(h_0 + h_1) & h_1 & 0 & \cdots & \\ 0 & h_1 & 2(h_1 + h_2) & h_2 & & \\ \vdots & & \ddots & \ddots & \ddots & \\ & & h_{n-3} & 2(h_{n-3} + h_{n-2}) & h_{n-2} & \\ & & & h_{n-2} & 2(h_{n-2} + h_{n-1}) & h_{n-1} \\ & & & & h_{n-1} & 2h_{n-1} \\ \end{bmatrix}}_{\boldsymbol{H}_{(n+1)\times(n+1)}} \underbrace{\begin{bmatrix} c_0 \\ c_1 \\ c_2 \\ \vdots \\ c_{n-2} \\ c_{n-1} \\ c_n \end{bmatrix}}_{\boldsymbol{C}_{(n+1) \times 1}} = \underbrace{\begin{bmatrix} 3(a_1 - a_0)/h_0 - 3f'(x_0) \\ 3(a_2 - a_1)/h_1 - 3(a_1 - a_0)/h_0 \\ 3(a_3 - a_2)/h_2 - 3(a_2 - a_1)/h_1 \\ \vdots \\ 3(a_{n-1} - a_{n-2})/h_{n-2} - 3(a_{n-2} - a_{n-3})/h_{n-3} \\ 3(a_n - a_{n-1})/h_{n-1} - 3(a_{n-1} - a_{n-2})/h_{n-2} \\ 3f'(x_n) - 3(a_n - a_{n-1})/h_{n-1} \end{bmatrix}}_{\boldsymbol{\beta}_{(n+1) \times 1}} \end{equation} $$

得到的 \(\boldsymbol{H}\) 还是一个 \((n+1) \times (n+1)\) 的方阵，可以解出 \(c_i, i \in \{0, \cdots, n-1\}\)。我们重载了 CubicSpline 的构造函数， 并调用成员函数 InitClampedSpline 来完成固定三次样条的构造。

4. 完