第二章 样条插值

分段低次插值具有收敛性好、数值稳定性高的优点，但光滑性较差。样条插值是一种既能保持分段低次插值的稳定性，又能获得较高光滑性的插值方法。

样条（Spline）一词来源于工程绘图中的样条曲线——一种具有弹性的细长木条，用以画出光滑曲线。数学上的样条函数正是模拟这种物理样条的特性而构造的。

定义 2.1（样条函数） 给定区间 $[a, b]$ 的一个分划 $\Delta: a = x_0 < x_1 < \cdots < x_n = b$，若函数 $S(x)$ 满足：

1. $S(x)$ 在每个小区间 $[x_i, x_{i+1}]$（$i = 0, 1, \ldots, n-1$）上是次数不超过 $m$ 的多项式； 2. $S(x)$ 在 $[a, b]$ 上具有 $m-1$ 阶连续导数。

则称 $S(x)$ 为关于分划 $\Delta$ 的 $m$ 次样条函数，$x_0, x_1, \ldots, x_n$ 称为样条节点。

$m$ 次样条函数空间记为 $S_m(\Delta)$，这是一个 $n+m$ 维的线性空间。

定义 2.2（三次样条插值） 给定函数 $f(x)$ 在节点 $x_0 < x_1 < \cdots < x_n$ 处的函数值 $y_i = f(x_i)$（$i = 0, 1, \ldots, n$），求三次样条函数 $S(x) \in S_3(\Delta)$ 满足： $$S(x_i) = y_i, \quad i = 0, 1, \ldots, n$$

在每个小区间 $[x_i, x_{i+1}]$ 上，$S(x)$ 是三次多项式，共 $4n$ 个待定系数。

约束条件：

1. 插值条件：$S(x_i) = y_i$，$n+1$ 个条件
2. 连续性：$S(x_i^-) = S(x_i^+)$，$n-1$ 个条件
3. 一阶导数连续：$S'(x_i^-) = S'(x_i^+)$，$n-1$ 个条件
4. 二阶导数连续：$S(x_i^-) = S(x_i^+)$，$n-1$ 个条件

共 $4n-2$ 个条件，还差 2 个条件，需要补充边界条件。

第一类边界条件（固定斜率）： $$S'(x_0) = f'_0, \quad S'(x_n) = f'_n$$

第二类边界条件（固定曲率）： $$S(x_0) = f_0, \quad S(x_n) = f_n$$

特别地，当 $f_0 = f_n = 0$ 时，称为自然边界条件，相应的样条称为自然样条。

第三类边界条件（周期条件）：当 $f(x)$ 是周期函数时 $$S'(x_0) = S'(x_n), \quad S(x_0) = S(x_n)$$

设 $S(x_i) = M_i$（称为弯矩），在每个小区间 $[x_i, x_{i+1}]$ 上，$S(x)$ 是线性函数： $$S(x) = M_i \frac{x_{i+1} - x}{h_i} + M_{i+1} \frac{x - x_i}{h_i}, \quad x \in [x_i, x_{i+1}]$$ 其中 $h_i = x_{i+1} - x_i$。 积分两次并利用插值条件，得到： $$S(x) = M_i \frac{(x_{i+1}-x)^3}{6h_i} + M_{i+1} \frac{(x-x_i)^3}{6h_i} + \left(y_i - \frac{M_i h_i^2}{6}\right)\frac{x_{i+1}-x}{h_i} + \left(y_{i+1} - \frac{M_{i+1} h_i^2}{6}\right)\frac{x-x_i}{h_i}$$ 由 $S'(x)$ 在内部节点连续，得到三弯矩方程： $$\mu_i M_{i-1} + 2M_i + \lambda_i M_{i+1} = d_i, \quad i = 1, 2, \ldots, n-1$$ 其中： - $\mu_i = \frac{h_{i-1}}{h_{i-1}+h_i}$ - $\lambda_i = \frac{h_i}{h_{i-1}+h_i} = 1 - \mu_i$ - $d_i = \frac{6}{h_{i-1}+h_i}\left(\frac{y_{i+1}-y_i}{h_i} - \frac{y_i-y_{i-1}}{h_{i-1}}\right) = 6f[x_{i-1}, x_i, x_{i+1}]$ 等距节点情形：$h_i = h$，则 $\mu_i = \lambda_i = \frac{1}{2}$，方程简化为： $$M_{i-1} + 4M_i + M_{i+1} = \frac{6}{h^2}(y_{i-1} - 2y_i + y_{i+1})$$ ==== 2.3.4 边界条件的处理 ==== 第一类边界条件： 由 $S'(x_0) = f'_0$ 和 $S'(x_n) = f'_n$，补充两个方程： $$2M_0 + M_1 = \frac{6}{h_0}\left(\frac{y_1-y_0}{h_0} - f'_0\right) = d_0$$ $$M_{n-1} + 2M_n = \frac{6}{h_{n-1}}\left(f'_n - \frac{y_n-y_{n-1}}{h_{n-1}}\right) = d_n$$ 得到三对角方程组： $$\begin{bmatrix} 2 & \lambda_0 & & &
\mu_1 & 2 & \lambda_1 & &
& \ddots & \ddots & \ddots &
& & \mu_{n-1} & 2 & \lambda_{n-1}
& & & \mu_n & 2 \end{bmatrix} \begin{bmatrix} M_0
M_1
\vdots
M_{n-1}
M_n \end{bmatrix} = \begin{bmatrix} d_0
d_1
\vdots
d_{n-1}
d_n \end{bmatrix}$$ 其中 $\lambda_0 = \mu_n = 1$。 第二类边界条件： 直接给定 $M_0 = f_0$，$M_n = f_n$，求解中间 $n-1$ 个未知数。 追赶法求解： 上述三对角方程组可用追赶法（Thomas算法）高效求解，计算复杂度为 $O(n)$。 ==== 2.3.5 三转角方程 ==== 另一种方法是设 $S'(x_i) = m_i$（称为转角），可导出三转角方程： $$\lambda_i m_{i-1} + 2m_i + \mu_i m_{i+1} = g_i, \quad i = 1, 2, \ldots, n-1$$ 其中： $$g_i = 3\left(\lambda_i \frac{y_i - y_{i-1}}{h_{i-1}} + \mu_i \frac{y_{i+1} - y_i}{h_i}\right)$$ ===== 2.4 B样条 ===== ==== 2.4.1 截断幂函数 ==== 定义 2.3（截断幂函数） $$x_+^m = \begin{cases} x^m, & x \geq 0
0, & x < 0 \end{cases}$$ 截断幂函数 $x_+^m$ 具有 $m-1$ 阶连续导数。 ==== 2.4.2 B样条的定义 ==== 定义 2.4（B样条） 给定节点 $t_0 \leq t_1 \leq \cdots \leq t_{m+1}$，$m$ 次 B样条 定义为： $$B_{i,m}(x) = (t_{i+m+1} - t_i)[t_i, t_{i+1}, \ldots, t_{i+m+1}](t - x)_+^m$$ 其中 $[t_i, \ldots, t_{i+m+1}](t-x)_+^m$ 表示对变量 $t$ 的 $m+1$ 阶差商。 递推公式（de Boor-Cox公式）： $$B_{i,0}(x) = \begin{cases} 1, & t_i \leq x < t_{i+1}
0, & \text{其他} \end{cases}$$ $$B_{i,m}(x) = \frac{x - t_i}{t_{i+m} - t_i} B_{i,m-1}(x) + \frac{t_{i+m+1} - x}{t_{i+m+1} - t_{i+1}} B_{i+1,m-1}(x)$$ ==== 2.4.3 B样条的性质 ==== 1. 局部支撑性：$B_{i,m}(x) = 0$ 当 $x \notin [t_i, t_{i+m+1}]$ 2. 非负性：$B_{i,m}(x) \geq 0$ 对所有 $x$ 3. 单位分解：$\sum_i B_{i,m}(x) = 1$ 4. 线性无关性：$B_{i,m}(x)$ 在其支撑集上线性无关 ==== 2.4.4 三次B样条 ==== 对于等距节点，三次B样条（$m=3$）的显式表达式为： $$B_{i,3}(x) = \frac{1}{6h^3} \begin{cases} (x-t_{i-2})^3, & x \in [t_{i-2}, t_{i-1})
h^3 + 3h^2(x-t_{i-1}) + 3h(x-t_{i-1})^2 - 3(x-t_{i-1})^3, & x \in [t_{i-1}, t_i)
h^3 + 3h^2(t_{i+1}-x) + 3h(t_{i+1}-x)^2 - 3(t_{i+1}-x)^3, & x \in [t_i, t_{i+1})
(t_{i+2}-x)^3, & x \in [t_{i+1}, t_{i+2})
0, & \text{其他} \end{cases}$$ 三次B样条插值： $$S(x) = \sum_{i=-1}^{n+1} c_i B_{i,3}(x)$$ ===== 2.5 样条逼近 ===== ==== 2.5.1 最小二乘样条逼近 ==== 给定数据点 $(x_i, y_i)$（$i = 0, 1, \ldots, N$），$N > n$，求 $m$ 次样条函数： $$S(x) = \sum_{j=-m}^{n-1} c_j B_{j,m}(x)$$ 使得： $$\Phi = \sum_{i=0}^{N} \omega_i [y_i - S(x_i)]^2$$ 最小。这导出一个关于系数 $c_j$ 的线性方程组。 ==== 2.5.2 光顺样条 ===== 考虑泛函： $$J(S) = \sum_{i=0}^{n} \omega_i [y_i - S(x_i)]^2 + \lambda \int_a^b [S(x)]^2 dx$$

其中 $\lambda > 0$ 是光滑参数：

1. $\lambda \to 0$：插值样条
2. $\lambda \to \infty$：线性函数（最小二乘直线）

最小化 $J(S)$ 得到的样条称为光顺样条或光滑样条（Smoothing Spline）。

定理 2.1（误差估计） 设 $f(x) \in C^4[a, b]$，$S(x)$ 是 $f(x)$ 的三次样条插值函数（满足第一类或第二类边界条件），则： $$\|f^{(k)} - S^{(k)}\|_\infty \leq C_k h^{4-k} \|f^{(4)}\|_\infty, \quad k = 0, 1, 2$$

其中 $h = \max h_i$，$C_0 = \frac{5}{384}$，$C_1 = \frac{1}{24}$，$C_2 = \frac{3}{8}$。

推论：当 $h \to 0$ 时，$S(x)$，$S'(x)$，$S(x)$ 分别一致收敛于 $f(x)$，$f'(x)$，$f(x)$。

例题 2.1 已知数据：

求满足自然边界条件的三次样条插值函数。

解：等距节点，$h = 1$，$y_i = 0$。

由自然边界条件：$M_0 = M_3 = 0$。

三弯矩方程（$i=1,2$）： $$M_0 + 4M_1 + M_2 = 6(y_0 - 2y_1 + y_2) = 0$$ $$M_1 + 4M_2 + M_3 = 6(y_1 - 2y_2 + y_3) = 0$$

即： $$4M_1 + M_2 = 0$$ $$M_1 + 4M_2 = 0$$

解得：$M_1 = M_2 = 0$。

因此 $S(x) = 0$（零函数）。这是显然的，因为所有插值点都是0。

例题 2.2 已知数据：

边界条件：$S'(0) = 1$，$S'(3) = 0$。求三次样条插值函数。

解：$h = 1$，先计算 $d_i$：

$d_0 = 6(f[0,1] - 1) = 6(2-1) = 6$

$d_1 = 6(f[0,1,2]) = 6 \cdot \frac{2-3-2}{2} = 6 \cdot (-1.5) = -9$

$d_2 = 6(f[1,2,3]) = 6 \cdot \frac{4-2-(-1)}{2} = 6 \cdot 1.5 = 9$

$d_3 = 6(0 - f[2,3]) = 6(0 - 2) = -12$

三对角方程组： $$\begin{cases} 2M_0 + M_1 = 6
\frac{1}{2}M_0 + 2M_1 + \frac{1}{2}M_2 = -9
\frac{1}{2}M_1 + 2M_2 + \frac{1}{2}M_3 = 9
M_2 + 2M_3 = -12 \end{cases}$$

解得：$M_0 = 6$，$M_1 = -6$，$M_2 = 6$，$M_3 = -9$。

在各小区间上的表达式可按公式写出。

基础练习

1. 什么是样条函数？三次样条插值需要满足哪些条件？

2. 证明三次样条函数空间 $S_3(\Delta)$ 的维数为 $n+3$。

3. 已知 $f(0) = 0$，$f(1) = 1$，$f(2) = 4$，求满足 $S(0) = S(2) = 0$ 的三次样条插值函数。

进阶练习

4. 证明B样条的单位分解性质：$\sum_i B_{i,m}(x) = 1$。

5. 设 $S(x)$ 是 $f(x)$ 的三次样条插值函数，证明：

 $$\int_a^b [S''(x)]^2 dx \leq \int_a^b [f''(x)]^2 dx$$
 等号当且仅当 $f(x) = S(x)$ 时成立。

6. 设 $f(x) \in C^2[a, b]$，$S_I(x)$ 是第一类边界条件下的样条插值，$S_{II}(x)$ 是自然样条插值。证明：

 $$\int_a^b [f''(x) - S''_{II}(x)]^2 dx = \int_a^b [f''(x)]^2 dx - \int_a^b [S''_{II}(x)]^2 dx + 2f''(b)S''_{II}(b) - 2f''(a)S''_{II}(a)$$

编程实践

7. 编写程序实现三弯矩法求解三次样条插值。

8. 用Python的SciPy库（scipy.interpolate.CubicSpline）实现三次样条插值，并与自己编写的程序进行对比。

1. 李庆扬等. 数值分析（第5版）. 清华大学出版社, 2008. 2. de Boor C. A Practical Guide to Splines. Springer, 2001. 3. Schumaker L L. Spline Functions: Basic Theory. Cambridge University Press, 2007.