张叶安的小站

本页面只读。您可以查看源文件，但不能更改它。如果您觉得这是系统错误，请联系管理员。
====== 第二章 样条插值 ======

===== 2.1 引言 =====

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

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

===== 2.2 样条函数的基本概念 =====

**定义 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.3 三次样条插值 =====

==== 2.3.1 三次样条插值问题 ====

**定义 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$ 个待定系数。

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

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

==== 2.3.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)$$

==== 2.3.3 三弯矩方程 ====

设 $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$ 是光滑参数：
  - $\lambda \to 0$：插值样条
  - $\lambda \to \infty$：线性函数（最小二乘直线）

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

===== 2.6 收敛性与误差估计 =====

**定理 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.7 典型例题 =====

**例题 2.1** 已知数据：
| $x$ | 0 | 1 | 2 | 3 |
|-----|---|---|---|---|
| $y$ | 0 | 0 | 0 | 0 |

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

**解**：等距节点，$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** 已知数据：
| $x$ | 0 | 1 | 2 | 3 |
|-----|---|---|---|---|
| $f(x)$ | 1 | 3 | 2 | 4 |

边界条件：$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$。

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

===== 2.8 习题 =====

**基础练习**

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.