Part 2 有理分式内插

本节我们讨论利用有理分式进行内插。当函数在某个区间内的变化行为比较剧烈时，多项式内插会出现剧烈震荡的行为。这时候利用有理分式进行内插可能会更为合适一些。

考虑一个分式，

$$\Phi^{(m,n)}(x) = \frac{P_m(x)}{Q_n(x)} \ ,$$

其中分子和分母分别是 $m$ 阶和 $n$ 阶的多项式。我们下面将假定这两个多项式是互素的，也就是说并不存在一个公共的多项式因子。

下面我们直接构造这样的有理函数。从 $i = 0, 1, \cdots$ 开始，我们用下面的次序来递推：

$$\begin{array}{c|ccc|l} i & x_i & y_i & \text{others} & \\ \hline 0 & x_0 & y_0 & & \\ 1 & x_1 & y_1 & \phi(x_0, x_1) & \\ 2 & x_2 & y_2 & \phi(x_0, x_2) & \phi(x_0, x_1, x_2) \\ 3 & x_3 & y_3 & \phi(x_0, x_3) & \phi(x_0, x_1, x_3) \quad \phi(x_0, x_1, x_2, x_3) \\ \vdots & \vdots & \vdots & & \vdots \end{array}$$

其中的各个函数 $\phi$ 可以按照下式进行递推地定义：

$$\begin{aligned} \phi(x_i, x_j) &= \frac{x_i - x_j}{y_i - y_j} \ , \\ \phi(x_i, x_j, x_k) &= \frac{x_j - x_k}{\phi(x_i, x_j) - \phi(x_i, x_k)} \ , \\ &\vdots \\ \phi(x_i, \cdots, x_l, x_m, x_n) &= \frac{x_m - x_n}{\phi(x_i, \cdots, x_l, x_m) - \phi(x_i, \cdots, x_l, x_n)} \ . \end{aligned}$$

显然，为了获得具体的数而不是发散的结果，我们应当尽可能选择函数单调的一个区间进行内插的构造。由于我们总是假设各个 $x_i$ 是不相同的，因此如果函数是单调的，那么上述各个构造的差值比 $\phi$ 就不会发散。一旦获得了这些系数，我们可以构造一个有理分式：

$$\Phi^{(n,n)}(x) = \frac{P_n(x)}{Q_n(x)} \ ,$$

它是两个 $n$ 次多项式的比。我们要求它能够经过 $(2n+1)$ 个点 $(x_i, y_i), \ i = 0, 1, \cdots, 2n$。事实上这个结果可以写成一个连分数：

$$\Phi^{(n,n)}(x) = y_0 + \frac{x - x_0}{\phi(x_0, x_1) + \frac{x - x_1}{\phi(x_0, x_1, x_2) + \ddots + \frac{x - x_{2n-1}}{\phi(x_0, \cdots, x_{2n})}}} \ .$$

可以证明这个分式恰好通过给定的 $2n+1$ 个点，即满足 $\Phi^{(n,n)}(x_i) = y_i, \ i = 0, 1, \cdots, 2n$。

作为一个例子，我们考虑前面曾经提及的 Runge 函数：

$$f(x) = \frac{1}{1 + 25x^2} \ .$$

我们选择其单调的区间，比如 $[0, 1]$ 来进行内插。假定知道下列三个点的数值： $x^2 = 0, 1/25, 1$，相应的函数值分别为： $1, 1/2, 1/26$。按照上面的构造，按照 $x^2$ 的有理分式为，

$$f(x) = \Phi^{(1,1)}(x^2) = 1 + \frac{x^2}{\phi(0, 1/25) + \frac{x^2 - 1/25}{\phi(0, 1/25, 1)}}$$

其中的几个系数，

$$\begin{aligned} \phi(0, 1/25) &= \frac{0 - 1/25}{1 - 1/2} = -\frac{2}{25} \ , \\\\ \phi(0, 1) &= \frac{0 - 1}{1 - 1/26} = -\frac{25}{26} \ , \\\\ \phi(0, 1/25, 1) &= \frac{1/25 - 1}{-2/25 + 26/25} = -1 \ . \end{aligned}$$

将这些数值带入前式我们发现，我们竟然可以完全重构 $f(x)$。

另外一个例子，让我们考察 $f(x) = \cot(x)$ 的计算问题。

Remez 算法可以用来计算这种近似，而这种算法与所谓的 Chebyshev 多项式密切相关。

Part 3 样条函数内插

这一节中，我们来讨论利用样条函数 (spline function) 来进行内插。

考虑一个区间 $[a, b]$ 的一个 $n$ 段分割： $a = x_0 < x_1 \cdots x_{n-1} < x_n = b$。位于区间中间的那些点 $x_i, \ i = 1, \cdots, n-1$ 称为节点 (knots)。我们考虑在该区间上的一个所谓三次样条函数 (cubic spline function) $S_\Delta: [a, b] \to \mathbb{R}$，它满足 $S_\Delta \in C^2[a, b]$ 并且在每个子区间 $[x_i, x_{i+1}]$ 上均为一个三次多项式并且它们在每个节点处函数值以及一二阶导数值都连续 (平滑衔接)。如果给定的一组函数值记为 $Y = \{y_0, y_1, \cdots, y_n\}$。我们把这样的三次样条差值函数记为 $S_\Delta(Y; \cdot)$，它满足 $S_\Delta(Y; x_i) = y_i, \ i = 0, 1, \cdots, n$。

上述信息还没有完全确定函数 $S_\Delta(Y; \cdot)$ 的性质。可以证明还有两个自由度可以调节。这两个自由度一般通过在两个端点处的边条件来加以限制。边条件大致可以分为下列三种：第一种是 Dirichlet 边条件；第二种是 Neumann 边条件；第三种是 周期边条件。可以证明，上述三种方案的任何一种都唯一地确定了一个三次样条差值函数 $S_\Delta(Y; \cdot)$。样条函数被广泛的应用是因为它实际上并不像高阶的多项式内插那样会呈现出剧烈震荡的 Runge 现象。

我们现在考虑某个区间 $[a, b]$ 上的 $m > 0$ 阶以下的所有阶导数连续且其 $m$ 阶导数平方可积的函数构成的函数空间，将其记为 $\mathcal{K}^m[a, b]$。我们同时会用 $\mathcal{K}_p^m[a, b]$ 来标记所有的 $m$ 阶以下导数都满足周期性的函数，即 $f^{(k)}(a) = f^{(k)}(b)$， $k = 0, \cdots, m-1$。对于任意的函数 $f \in \mathcal{K}^2[a, b]$，我们可以定义一个函数的模：

$$||f|| = \int_a^b dx |f''(x)|^2 \ . \tag{b}$$

我们知道一个函数 $f(x)$ 的曲率 $R(x) = f''(x)(1+f'(x)^2)^{-3/2}$。如果在一个区间上 $f'(x)^2$ 比起 $1$ 来说要小很多的话，那么函数的曲率几乎就等于 $f''(x)$。因此，上面这个模从某种意义上是衡量了函数在一个区间上曲率模方的大小。可以证明，三次样条函数实际上是使得这个模最小的函数。具体来说，我们不加证明地引用下列定理，它一般称为 最小模定理。从直观上说，它是“最光滑的”函数。

定理 3.1 实区间 $[a, b]$ 的任一 $n$ 段分割 $\Delta: a = x_0 < x_1 \cdots x_{n-1} < x_n = b$ 上的三次样条函数 $S_\Delta(Y; x_i) = y_i$， $i = 0, 1, \cdots, n$。对任给的 $f \in \mathcal{K}^2[a, b]$ 且 $f(x_i) = y_i$，我们一定有 $||f|| \ge ||S_\Delta(Y; \cdot)||$。其中的函数模由式 (b) 所定义。事实上可以证明，

$$||f - S_\Delta(Y; \cdot)||^2 = ||f||^2 - ||S_\Delta(Y; \cdot)||^2 \ge 0 \ ,$$

只要函数 $f$ 满足下列条件之一即可：

1. $S_\Delta''(Y; a) = S_\Delta''(Y; b) = 0$;
2. $f \in \mathcal{K}_p^2[a, b]$，$S_\Delta(Y; \cdot)$ 是周期的;
3. $f'(a) = S_\Delta'(Y; a)$， $f'(b) = S_\Delta'(Y; b)$.

并且上面的任何一个条件都唯一地确定了样条函数 $S_\Delta(Y; \cdot)$。

下面我们具体讨论如何计算确定一个区间上的三次样条函数。我们会发现，它实际上对应于求解一个三对角矩阵的线性方程组。

首先，对于任意的一个 $n$ 段分割 $\Delta: a = x_0 < x_1 \cdots x_{n-1} < x_n = b$，我们首先引入下列记号：

$$h_{j+1} = x_{j+1} - x_j \ , j = 1, \cdots, n \ ,$$

它标志了每一段的长度。其次，我们注意到对于三次样条函数 $S_\Delta(Y; \cdot)$ 来说，它的二次导数 $S_\Delta''(Y; \cdot)$ 在每一段上面必定是一个一次函数。我们将其在各个节点处的数值称为三次样条函数的矩 (moments) 并记为，

$$M_j = S_\Delta''(Y; x_j) \ , j = 0, 1, \cdots, n-1, n \ .$$

对于每一小段 $x \in [x_j, x_{j+1}]$， $j = 0, \cdots, n-1$ 上的线性函数 $S_\Delta''(Y; \cdot)$，我们可以验明，

$$S_\Delta''(Y; x) = M_j \frac{x_{j+1} - x}{h_{j+1}} + M_{j+1} \frac{x - x_j}{h_{j+1}} \ .$$

以此式为出发点我们可以积分一次获得，

$$S_\Delta'(Y; x) = -M_j \frac{(x_{j+1} - x)^2}{2h_{j+1}} + M_{j+1} \frac{(x - x_j)^2}{2h_{j+1}} + A_j \ , \tag{c}$$

进一步再积分一次得到，

$$S_\Delta(Y; x) = M_j \frac{(x_{j+1} - x)^3}{6h_{j+1}} + M_{j+1} \frac{(x - x_j)^3}{6h_{j+1}} + A_j(x - x_j) + B_j \ . \tag{d}$$

利用条件 $S_\Delta(Y; x_j) = y_j$， $S_\Delta(Y; x_{j+1}) = y_{j+1}$ 我们就可以求出积分常数 $A_j$ 和 $B_j$，其结果如下：

$$\begin{cases} B_j = y_j - M_j \frac{h_{j+1}^2}{6} \ , \\ A_j = \frac{y_{j+1} - y_j}{h_{j+1}} - \frac{h_{j+1}}{6}(M_{j+1} - M_j) \ . \end{cases} \tag{e}$$

公式 (d) 和公式 (e) 告诉我们，一个任意 $n$ 点分割 $\Delta$ 上的三次样条函数 $S_\Delta(Y; \cdot)$ 完全由其在各个节点 $x_j$ 处的矩 $M_j$ 所完全确定。于是，剩下的任务就是根据已知的条件确定样条函数在各节点处的矩。

我们现在还没有利用 $S_\Delta'(Y; x_j)$ 在各个节点处平滑衔接的条件。将公式 (e) 带入公式 (d) 并要求 $x_j$ 两侧的导数平滑衔接 $S_\Delta'(Y; x_j + 0) = S_\Delta'(Y; x_j - 0)$，我们得到，

$$\frac{h_j}{6} M_{j-1} + \frac{h_j + h_{j+1}}{3} M_j + \frac{h_{j+1}}{6} M_{j+1} = \frac{y_{j+1} - y_j}{h_{j+1}} - \frac{y_j - y_{j-1}}{h_j} \ . \tag{f}$$

这个公式中的 $j = 1, 2, \cdots, n-1$ 给出了 $(n-1)$ 个方程。但是我们的未知数 $M_j$ 共有 $(n+1)$ 个。我们仍然需要两个条件才可能完全确定所有的矩 $M_j$ 进而确定三次样条函数。这两个额外的条件恰恰由我们定理 3.1 中所提及的，函数在端点处的第一类、第二类或周期性条件给出。

公式 (f) 已经体现出典型的三对角矩阵所对应的线性方程组的形式。我们完全可以利用第二章第 7 节中介绍的方法进行数值求解。

更新日志

2026/10/7 03:54

查看所有更新日志

• 98079-修正于 2026/10/7
• 76ef6-computational-physics于 2026/9/22

版权所有

版权归属：nicostore-mathematica

许可证：署名-非商业性-相同方式共享 4.0 国际 (CC-BY-NC-SA-4.0)

贡献者: nicostore-mathematica