4. 有限差分系数设计
在数值微分中，差分格式的系数并不是随意凑出来的，而是由一套严格的代数方法确定。给定一个节点模板（例如三点、五点），我们希望用这些节点上的函数值线性组合去逼近某阶导数，而组合系数必须满足一组由泰勒展开或多项式插值导出的线性方程。本节围绕这套系数求解方法展开，先讲清原理，再给出可复现的推导与计算过程，最后落到二阶至四阶常用格式的系数表与使用判断。
一、系数设计的本质：让误差项尽可能高阶地消失
设我们要在节点 x0 处逼近 f′(x0)，取等距步长 h，模板节点为 x0+jh（j 取一组整数）。构造线性组合
D(f) = (1/h) · Σj aj f(x0+jh)
把每个 f(x0+jh) 在 x0 处做泰勒展开：
f(x0+jh) = f0 + jh·f0' + (j²h²/2)·f0'' + (j³h³/6)·f0''' + (j⁴h⁴/24)·f0⁗ + …代入线性组合并按 h 的幂次归并，得到
D(f) = (1/h)·[ (Σ a_j) f0 + h(Σ j a_j) f0' + (h²/2)(Σ j² a_j) f0''
+ (h³/6)(Σ j³ a_j) f0''' + … ]要让 D(f) 逼近 f0′，就必须让 f0 的系数为零、f0′ 的系数为 1，并尽可能让更高阶导数的系数为零。于是得到矩条件：
Σ a_j = 0
Σ j·a_j = 1
Σ j²·a_j = 0
Σ j³·a_j = 0
…模板中有 m 个节点就有 m 个待定系数，可以满足前 m 个矩条件。前两个条件保证一阶精度，之后每多满足一个条件，精度提高一阶。这就是“系数设计”的全部代数骨架：节点位置决定矩阵，矩条件决定右端项，解线性方程组得到系数。若目标是二阶导数，则把右端项改为 Σ j·aj=0、Σ j²aj=2，其余条件不变。
二、从泰勒展开到线性方程组：三点格式的完整推导
取对称模板 j = −1, 0, 1，待定系数 a−1、a0、a1。三个未知数对应三个矩条件：
a₋₁ + a₀ + a₁ = 0
−a₋₁ + a₁ = 1
a₋₁ + a₁ = 0由第二、三式直接得 a1 = 1/2，a−1 = −1/2，再代入第一式得 a0 = 0。于是
f'(x0) ≈ [f(x0+h) − f(x0−h)] / (2h)这就是中心差分。由于模板对称，所有偶数阶矩条件自动满足（Σ j2kaj=0 对奇数阶导数目标成立），因此误差项中 h² 项消失，截断误差为 O(h²)。这一点很关键：对称模板天然消去偶数阶误差项，是中心格式精度高的根本原因，而不是巧合。
若改用前向模板 j = 0, 1, 2，同样三个条件解得 a0 = −3/2，a1 = 2，a2 = −1/2，即
f'(x0) ≈ [−3f(x0) + 4f(x0+h) − f(x0+2h)] / (2h)这是二阶精度的前向格式。它只用右侧节点，适合边界处无法取到左侧点的场景，但误差常数通常大于同阶中心格式。
三、多项式插值视角：同一系数的另一种来源
系数也可以由插值多项式求导得到。取模板上 m 个节点构造拉格朗日插值多项式 L(x)，则 f′(x0) ≈ L′(x0)。对等距节点，L′(x0) 展开后每个节点值的权重恰好就是上面解出的 aj/h。两种方法殊途同归：泰勒法直接给矩条件，插值法给出权重的显式表达。实践中，泰勒法更适合手推和验证，插值法更适合程序化生成任意阶格式。例如用插值法可写出五点中心一阶导格式的权重为 (1, −8, 0, 8, −1)/12，对应
f'(x0) ≈ [f(x0−2h) − 8f(x0−h) + 8f(x0+h) − f(x0+2h)] / (12h)其截断误差为 O(h⁴)。注意中间节点 x0 的系数为 0，这是对称模板的必然结果，也意味着该格式不消耗中心点的函数值。
四、二阶至四阶常用格式的系数与误差对照
下表汇总了工程中最常用的几组格式。表中“精度阶”指截断误差的主项阶数，系数已按 1/h 归一化，使用时直接乘以对应函数值再除以 h 即可。
| 格式 | 节点模板（相对 x₀） | 系数 aj | 精度阶 | 典型用途 |
|---|---|---|---|---|
| 一阶导·前向 | 0, 1 | −1, 1 | O(h) | 快速估算、边界单侧 |
| 一阶导·中心 | −1, 0, 1 | −1/2, 0, 1/2 | O(h²) | 内点通用 |
| 一阶导·前向二阶 | 0, 1, 2 | −3/2, 2, −1/2 | O(h²) | 左边界 |
| 一阶导·五点中心 | −2, −1, 0, 1, 2 | 1/12, −2/3, 0, 2/3, −1/12 | O(h⁴) | 高精度内点 |
| 二阶导·中心 | −1, 0, 1 | 1, −2, 1 | O(h²) | 曲率、扩散项 |
| 二阶导·五点中心 | −2, −1, 0, 1, 2 | −1/12, 4/3, −5/2, 4/3, −1/12 | O(h⁴) | 高精度二阶导 |
需要强调，精度阶高不等于任何步长下都更准。截断误差是 O(hp) 的渐近描述，当 h 不够小、函数高阶导数较大时，高阶格式的误差常数可能反而更大；同时高阶格式模板更宽，在边界附近需要更多单侧处理。因此选择格式时要同时看精度阶、模板宽度和边界可用性。
五、实际计算：用程序化方式生成任意格式系数
手推适合理解，但工程中更稳妥的做法是用程序解矩条件方程组，避免符号错误。下面给出一个可直接运行的 Python 片段，输入节点偏移列表和目标导数阶，输出归一化系数与精度阶。它把“节点位置决定矩阵、矩条件决定右端项”这一逻辑直接编码。
import numpy as np
def fd_coeffs(offsets, deriv=1):
"""offsets: 节点相对偏移，如 [-1,0,1]；deriv: 目标导数阶"""
m = len(offsets)
A = np.array([[j**k for j in offsets] for k in range(m)], dtype=float)
b = np.zeros(m)
b[deriv] = np.math.factorial(deriv) # 目标阶导数的矩条件右端项
a = np.linalg.solve(A, b)
# 精度阶：第一个未被满足的矩条件阶数减 1
order = m - deriv
return a, order
print(fd_coeffs([-1, 0, 1], 1)) # 中心差分
print(fd_coeffs([0, 1, 2], 1)) # 前向二阶
print(fd_coeffs([-2,-1,0,1,2], 1))# 五点中心以三点中心模板为例，矩阵 A 的三行分别是 j0、j1、j2 在 j=−1,0,1 处的取值，右端向量为 (0,1,0)，解出 (−0.5, 0, 0.5)，与手推一致。对五点中心一阶导，解出 (1/12, −2/3, 0, 2/3, −1/12)，精度阶为 4。运行后应核对两点：系数之和是否为零（保证常数项消去），以及 Σ j·aj 是否等于目标阶乘（保证主项系数正确）。这两条是最快的自检手段。
六、边界处理与格式选择的实践判断
内点用中心格式最自然，但区间端点无法取到两侧节点，必须改用单侧模板。常见做法是：左端点用前向二阶格式，右端点用其镜像的后向格式，内点用中心格式。这样整体仍保持二阶精度，代价是边界误差常数偏大。若追求四阶整体精度，则边界也需配套的四阶单侧格式，例如左端五点前向一阶导系数为 (−25/12, 4, −3, 4/3, −1/4)，模板宽度增加到五个点，实现复杂度明显上升。
选择格式时可参考以下判断顺序：先看精度需求与步长量级，若 h 已很小且函数光滑，二阶中心格式往往足够；再看边界可用节点数，节点不足时降级为单侧低阶格式；最后看函数高阶导数的量级，若 f‴ 很大，二阶格式的误差常数可能主导，此时提升到四阶格式更划算。需要提醒的是，高阶格式对噪声更敏感——差分本身是放大高频误差的运算，模板越宽、系数绝对值越大，噪声放大越明显。这正是后续噪声抑制要处理的问题，但在系数设计阶段就应意识到：精度阶与抗噪性之间存在权衡，不能只看阶数。
把本节方法落到一句话：给定节点模板，用泰勒展开写出矩条件，解一个范德蒙德型线性方程组，就得到该模板下唯一的一组差分系数；对称模板自动获得更高精度，单侧模板用于边界，而格式的最终选择要综合精度、边界条件与噪声敏感性三方面判断。