引言
数值计算的几乎每个环节都是”用简单函数代替复杂函数”:数值微积分把被积函数换成插值多项式,有限元把解空间换成分片多项式空间,嵌入式控制把昂贵的传感器模型换成查表加插值。逼近论为这一切提供误差语言与最优性判据。本文围绕三个基本问题展开:给定节点上的函数值,如何插值?光滑性要求高时如何分段拼接?在全体 $n$ 次多项式中,谁离目标最近?三者的答案——Lagrange 插值、三次样条与 Chebyshev 最佳一致逼近——之间由 Runge 现象与 Chebyshev 点的故事紧密串联。
理论基础
插值与误差公式
给定互异节点 $x_0, \dots, x_n$ 上的函数值 $f_i$,存在唯一 $n$ 次插值多项式 $p_n$。若 $f \in C^{n+1}[a,b]$,点态误差为
$$f(x) - p_n(x) = \frac{f^{(n+1)}(\xi_x)}{(n+1)!}\prod_{i=0}^{n}(x - x_i), \qquad \xi_x \in (a, b)$$
误差由两部分相乘:高阶导数的大小(函数本身的光滑性)与节点多项式 $\omega_{n+1}(x) = \prod (x-x_i)$ 的振幅(节点的选择)。等距节点下 $\omega_{n+1}$ 在区间端点附近急剧放大,这就是著名的 Runge 现象:对 $f(x) = 1/(1+25x^2)$ 在 $[-1,1]$ 上等距加密插值,端部误差不但不收敛反而指数发散。更精确地说,插值算子的 Lebesgue 常数 $\Lambda_n$(误差放大因子的上界)在等距节点下随 $n$ 指数增长。
Chebyshev 点与最佳逼近
Chebyshev 节点把节点多项式的最大振幅压到最小:取
$$x_j = \cos\frac{(2j+1)\pi}{2n+2}, \qquad j = 0, 1, \dots, n$$
即 Chebyshev 多项式 $T_{n+1}$ 的零点,则 $\max_{x\in[-1,1]}|\omega_{n+1}(x)| = 2^{-n}$ 为所有首一 $n+1$ 次多项式中的最小值,Lebesgue 常数仅以对数速度增长。Chebyshev 多项式由三项递推生成:
$$T_{k+1}(x) = 2x,T_k(x) - T_{k-1}(x), \qquad T_0 = 1,\quad T_1 = x$$
进一步,在全体 $n$ 次多项式中使 $|f - p|_\infty$ 最小者称为最佳一致逼近多项式 $p_n^*$,Chebyshev 交错定理给出充要判据:$p = p_n^*$ 当且仅当误差 $f – p$ 在至少 $n+2$ 个点上交错地达到 $\pm|f-p|_\infty$。该定理是 Remez 交换算法的理论基础,数学函数库中 $\sin$、$\exp$ 的硬件级实现即用其产生的极小极大多项式。
样条:用分段换稳定
另一条思路是放弃单一高次多项式,改用分段低次多项式并控制拼接光滑性。三次样条在每个子区间上取三次多项式,整体 $C^2$ 连续。Holladay 定理指出:在所有插值同一组数据的 $C^2$ 函数中,自然三次样条使曲率能量 $\int (f'')^2,\mathrm{d}x$ 最小——物理上对应柔性木条穿过固定钉点的平衡形态,”样条”一词即源于绘图工的木条。其构造归结为求解一个对称三对角线性方程组(未知量为节点二阶导数 $M_i$),代价 $O(n)$,逼近阶 $O(h^4)$。B-样条基由 de Boor–Cox 递推定义,具有局部支集与单位分解性,是 CAD 中 NURBS 曲线的底层积木。
数值方法与实现
import numpy as np
def bary_weights(x):
"""重心插值权重"""
n = len(x)
w = np.ones(n)
for j in range(n):
w[j] = 1.0/np.prod(x[j] - np.delete(x, j))
return w
def bary_eval(x, w, f, xx):
"""重心形式求值(第二形式)"""
num = np.zeros_like(xx); den = np.zeros_like(xx)
for j in range(len(x)):
diff = xx - x[j]
term = w[j]/diff
num += term*f[j]; den += term
num[diff == 0.0] = f[j]; den[diff == 0.0] = 1.0
return num/den
f = lambda x: 1.0/(1.0 + 25.0*x**2) # Runge 函数
xx = np.linspace(-1.0, 1.0, 2001)
for n in (10, 20, 40):
xe = np.linspace(-1.0, 1.0, n + 1) # 等距节点
xc = np.cos((2*np.arange(n + 1) + 1)*np.pi/(2*n + 2)) # Chebyshev 零点
for name, x in (('等距', xe), ('Chebyshev', xc)):
p = bary_eval(x, bary_weights(x), f(x), xx)
print(f"n={n:2d} {name:>9s} 最大误差 {np.max(np.abs(p - f(xx))):.3e}")
def natural_spline(x, y, xx):
"""自然三次样条:三对角方程求节点二阶导 M_i"""
n = len(x) - 1
h = np.diff(x)
A = np.zeros((n + 1, n + 1)); rhs = np.zeros(n + 1)
A[0, 0] = A[n, n] = 1.0
for i in range(1, n):
A[i, i-1] = h[i-1]; A[i, i] = 2.0*(h[i-1] + h[i]); A[i, i+1] = h[i]
rhs[i] = 6.0*((y[i+1] - y[i])/h[i] - (y[i] - y[i-1])/h[i-1])
M = np.linalg.solve(A, rhs)
idx = np.clip(np.searchsorted(x, xx) - 1, 0, n - 1)
hi = x[idx + 1] - x[idx]
a = (x[idx + 1] - xx)/hi; b = (xx - x[idx])/hi
return (a*y[idx] + b*y[idx + 1]
+ ((a**3 - a)*M[idx] + (b**3 - b)*M[idx + 1])*hi**2/6.0)
xs = np.linspace(-1.0, 1.0, 9)
err = np.max(np.abs(natural_spline(xs, f(xs), xx) - f(xx)))
print("自然样条(8 段)最大误差:", f"{err:.3e}")
运行可看到:等距插值误差随 $n$ 增大而爆炸,Chebyshev 插值几何收敛,而仅 8 段的自然样条误差已与高次 Chebyshev 插值相当——分段策略用极低的自由度代价买到了稳定性。
工程案例分析
嵌入式控制系统中的传感器标定是逼近论的隐形战场:热电偶、压力传感器的非线性响应曲线在出厂时离散标定,运行时依靠分段线性或三次样条插值实时换算。选型时单片机算力与存储受限,等距节点加分段低次多项式成为事实标准——这本质上是”用样条对抗 Runge”的工程直觉。CAD/CAM 领域则是 B-样条的天下:汽车覆盖件曲面用 NURBS 表示,求交、光顺、等距面计算全部建立在 de Boor 算法的局部求值上。
| 逼近方式 | 自由度来源 | 整体光滑性 | 收敛保障 | 典型用途 |
|---|---|---|---|---|
| 等距高次插值 | 次数 $n$ | $C^\infty$ | 无(Runge 发散) | 理论分析 |
| Chebyshev 插值 | 次数 $n$ | $C^\infty$ | 解析函数几何收敛 | 谱方法、函数逼近 |
| 三次样条 | 分段数 $N$ | $C^2$ | $O(h^4)$ | CAD、数据拟合 |
| 最佳一致逼近 | 次数 $n$ | $C^\infty$ | 交错定理最优 | 数学库函数实现 |
谱方法把 Chebyshev 插值推到偏微分方程求解:微分矩阵作用在 Chebyshev–Gauss–Lobatto 点上,对流体力学槽道流等光滑解问题,十几个模态即可达到有限差分数千网格的精度。信号处理中的 Parks–McClellan 算法是另一例经典应用:它把最佳一致逼近的 Remez 交换思想搬到 FIR 滤波器设计,使通带与阻带波纹等幅交错,按交错定理这就是给定阶数下的最优滤波器——逼近论的一条抽象定理,就这样住进了每一片 DSP 芯片。
结论与展望
逼近论的核心教训是:逼近的好坏由”函数的光滑性”与”节点/基函数的分布”共同决定,自由度多并不等于精度高。Chebyshev 思想(把节点推向端部、让误差等幅交错)与样条思想(局部支撑、以拼接换稳定)构成现代数值分析的两大支柱。当前活跃方向包括有理逼近(AAA 算法对极点与奇性的自适应捕捉)、深度网络作为逼近器的速率理论,以及插值与再生核 Hilbert 空间方法在不确定性量化中的融合。
发表回复