微分方程数值解

微分方程数值解用有限次计算近似微分方程的解。对于初值问题，它通常从已知状态出发，以小时间步不断推进，得到离散时刻的近似值。算法不仅要给出一条曲线，还需要说明步长改变时误差怎样变化，以及计算会不会制造原方程没有的振荡或增长。

已知现在的温差，怎样估计下一刻

考虑一个经过无量纲化的冷却问题： $y$ 是物体相对环境的温差比例，初值为 1；时间尺度取为使方程恰好成为 $y^{'} = - y, y (0) = 1.$ 它的精确解是 $y (t) = e^{- t}$ 。这里特意选择能求解析解的例子，是为了逐项检查算法；对于更复杂的方程，同样的推进步骤仍然可以使用。

由导数的含义，在当前时刻 $t_{n}$ 附近，可把函数近似为切线： $y (t_{n} + h) \approx y (t_{n}) + h f (t_{n}, y (t_{n})) .$ 其中 $y^{'} = f (t, y)$ ， $h > 0$ 是步长。实际计算不知道真值 $y (t_{n})$ ，就用已经算出的 $Y_{n}$ 代替，得到前向欧拉法： $t_{n + 1} = t_{n} + h, Y_{n + 1} = Y_{n} + h f (t_{n}, Y_{n}) .$ 大写 $Y_{n}$ 表示算法值，小写 $y (t_{n})$ 表示真解。Driscoll、Braun《Fundamentals of Numerical Computation》§6.2

取 $h = 0.5$ 。第一步用斜率 −1，从 $(0, 1)$ 推到 $Y_{1} = 1 + 0.5 (- 1) = 0.5$ 。第二步要重新计算斜率，此时算法认为状态为 0.5，因此斜率为 −0.5，得到 $Y_{2} = 0.5 + 0.5 (- 0.5) = 0.25$ 。而 $y (1) = e^{- 1} \approx 0.3678794412$ ，计算偏低约 0.117879。

上图用两段切线近似下降的指数曲线，下图用预测中点的斜率完成第一步，得到比欧拉更接近精确曲线的终点 — 金线是真解，青线连接数值状态；中点法先向前试半步，粉色点与短线表示试探状态及其斜率。

本例真解向上弯曲，起点切线落在真曲线下方，所以第一步低估并不意外。第二步又从已经偏低的位置继续出发，显示了误差会被后续步骤带着走。

用中途的斜率，为什么更准确

欧拉法用起点斜率代表整步变化。显式中点法先以欧拉方式预测半步的位置，再计算半步斜率： $k_{1} = f (t_{n}, Y_{n}), k_{2} = f (t_{n} + \frac{h}{2}, Y_{n} + \frac{h}{2} k_{1}), Y_{n + 1} = Y_{n} + h k_{2} .$ 这里 $k_{1}, k_{2}$ 都是变化率，还没有乘整步长度 $h$ 。

在第一步， $k_{1} = - 1$ ，预测中点状态是 $1 + 0.25 (- 1) = 0.75$ ，于是 $k_{2} = - 0.75$ ，终点为 $Y_{1} = 1 - 0.5 \times 0.75 = 0.625$ 。第二步从 0.625 出发，预测中点为 0.46875，得到 $Y_{2} = 0.625 - 0.5 \times 0.46875 = 0.390625$ 。在 $t = 1$ 的误差降为约 0.022746。

为什么中途试探有效？对本例直接把递推展开： $Y_{n + 1} = (1 - h + \frac{h^{2}}{2}) Y_{n} .$ 而精确推进一小步是 $y (t + h) = e^{- h} y (t)$ ，其中 $e^{- h} = 1 - h + \frac{h^{2}}{2} - \frac{h^{3}}{6} + \dots .$ 欧拉因子 $1 - h$ 只匹配到一次项，中点因子还匹配了二次项，单步遗漏便从二次量变成三次量。这给出了本例中精度改善的直接原因。

四次试探组成一个整步

经典四阶 Runge–Kutta 法，简称 RK4，在每一步使用四个斜率： $\begin{matrix} k_{1} & = f (t_{n}, Y_{n}), \\ k_{2} & = f (t_{n} + h / 2, Y_{n} + h k_{1} / 2), \\ k_{3} & = f (t_{n} + h / 2, Y_{n} + h k_{2} / 2), \\ k_{4} & = f (t_{n} + h, Y_{n} + h k_{3}), \\ Y_{n + 1} & = Y_{n} + \frac{h}{6} (k_{1} + 2 k_{2} + 2 k_{3} + k_{4}) . \end{matrix}$ 这四次试探共同确定一个整步，不是依次走了四个完整时间步。中间状态主要用于估计斜率，也不应当作已经确认的真解取样。Driscoll、Braun，§6.4.3，Definition 6.4.2

对 $h = 0.5, Y_{0} = 1$ ，四个斜率依次为 $k_{1} = - 1, k_{2} = - 0.75, k_{3} = - 0.8125, k_{4} = - 0.59375.$ 所以第一步为 $Y_{1} = 1 + \frac{0.5}{6} (- 1 - 1.5 - 1.625 - 0.59375) = 0.6067708333 \dots .$ 对方程 $y^{'} = - y$ ，这一算法的每步乘数恰为 $P_{4} (- h) = 1 - h + \frac{h^{2}}{2} - \frac{h^{3}}{6} + \frac{h^{4}}{24} .$ 第二步再乘同一个数，得到 $Y_{2} \approx 0.3681708442$ ，距 $e^{- 1}$ 约 0.000291403。虽然高阶法每步计算更多次斜率，但同一步长下可以显著减少误差。

单步误差为什么会累积

为了区分两种误差，假想每一步都从真解位置出发。欧拉法的单步缺陷为 $d_{n} = y (t_{n + 1}) - y (t_{n}) - h f (t_{n}, y (t_{n})) .$ 由 Taylor 公式，若这段真解满足 $| y^{″} (t) | \leq M$ ，则 $| d_{n} | \leq M h^{2} / 2$ 。有些教材把 $d_{n} / h$ 称为局部截断误差；这里使用未除以步长的缺陷，避免把两种阶数混在一起。

真实计算并不是每一步都回到真解。设全局误差 $e_{n} = Y_{n} - y (t_{n})$ ，并假设 $f$ 对状态满足 Lipschitz 界 $| f (t, a) - f (t, b) | \leq L | a - b |$ ，范围覆盖相关的真解和数值状态。两种推进式相减，得到 $e_{n + 1} = e_{n} + h (f (t_{n}, Y_{n}) - f (t_{n}, y (t_{n}))) - d_{n} .$ 取绝对值便有 $| e_{n + 1} | \leq (1 + h L) | e_{n} | + \frac{M h^{2}}{2} .$ 第一项传播已有误差，第二项加入新误差。从精确初值 $e_{0} = 0$ 出发，反复代入并求几何级数：当 $L > 0$ 时， $| e_{n} | \leq \frac{M h}{2 L} ((1 + h L)^{n} - 1) \leq \frac{M h}{2 L} (e^{L T} - 1), n h \leq T .$ 最后一步用了 $1 + h L \leq e^{h L}$ 。若 $L = 0$ ，直接求和得到 $| e_{n} | \leq M T h / 2$ 。因此在固定时间区间上，欧拉法的全局误差是一阶量 $O (h)$ ，尽管每次新增的缺陷是 $O (h^{2})$ 。

在足够光滑且满足相应稳定性条件的初值问题上，中点法和 RK4 的单步缺陷分别为 $O (h^{3})$ 、 $O (h^{5})$ ，全局误差分别为 $O (h^{2})$ 、 $O (h^{4})$ 。如果方程有不光滑点或事件造成的突变，需要在那些位置重新处理，不能无条件照用这个阶数。

把步长减半，是一项可以执行的检查

仍计算到 $t = 1$ ，下表列出与精确解比较的绝对误差。

步长 h	欧拉法	显式中点法	RK4
0.5	0.11787944	0.02274556	0.0002914030
0.25	0.05147319	0.00464959	0.0000147582
0.125	0.02427053	0.00105380	0.0000008308
0.0625	0.01180531	0.00025110	0.0000000493

当步长已经足够小时，再减半通常把三种误差分别缩小到约 $1 / 2, 1 / 4, 1 / 16$ 。大步长处的比例还会受更高阶项影响，因而不必精确等于这些数。

双对数坐标上欧拉、中点、RK4误差随步长减小下降，三条线的渐近斜率分别为一二四 — 横轴、纵轴都按对数刻度绘制。读者可将表中任意两组数取对数，核对线段的斜率。

没有解析解时，可比较 $h, h / 2, h / 4$ 的计算，并检查守恒量、非负性或已知极限情况。相邻两次结果接近是证据之一，但还需确认它们不是共同落在一个失稳或错误的算法分支上。自适应方法通常把局部误差与 $a t o l + r t o l | Y |$ 一类尺度比较，据此调节步长；这不等同于对全部时间点的全局误差给出同样大小的严格保证。

真解稳定，算法却可能振荡

现在只把衰减速率改为 $y^{'} = - λ y$ ，其中 $λ > 0$ 。真解是 $e^{- λ t}$ ，从正初值出发始终为正并趋于零。欧拉法则给出 $Y_{n + 1} = (1 - λ h) Y_{n} .$ 因此算法会衰减的条件是 $| 1 - λ h | < 1$ ，即 $0 < λ h < 2$ 。若还希望每一步保持非负，就要 $λ h \leq 1$ 。在 $1 < λ h < 2$ 时，幅值虽在衰减，符号却交替，产生原解没有的假振荡；在 $λ h > 2$ 时，假振荡还会放大。在边界 $λ h = 2$ ，幅值不再衰减，也没有模拟出真解的长期行为。

同一衰减方程用三种欧拉步长，乘数依次为0.5负0.5负1.5，表现为正值衰减、交替衰减、交替放大 — 三幅均为 λ=20。金线是真解，青线连接欧拉数值点；横轴是时间，点的疏密随步长改变。

后向欧拉法把斜率放在未知的下一时刻： $Y_{n + 1} = Y_{n} + h f (t_{n + 1}, Y_{n + 1})$ 。对本例可解出 $Y_{n + 1} = \frac{Y_{n}}{1 + λ h} .$ 任何正步长都给出正的衰减因子，但步长很大时仍可能不准确。对一般非线性方程，每一步还需求解一个代数方程；这里的方便除法是线性例子的特殊情况。Driscoll、Braun，Absolute stability

如果系统同时含有很快的衰减模式和缓慢变化，显式方法可能被快速模式限制步长，这称为刚性问题的一种典型表现。例如 $y^{'} = - 100 (y - \cos t) - \sin t, y (0) = 1$ 的精确解是缓慢变化的 $\cos t$ ；但相对于这条解的小扰动满足 $e^{'} = - 100 e$ 。前向欧拉的误差传播因子为 $1 - 100 h$ ，仍需 $h < 0.02$ 才衰减。只看目标曲线平缓，不能断定大步长就安全。

在SIR模型中，可检查总人数守恒与非负性；在弹簧振子模型中，可检查能量输入和损耗；在经典捕食者-猎物模型中，可检查守恒量是否漂移。这些检查针对计算本身。即使数值误差已经很小，假设机制和参数是否符合对象，仍是数学建模要回答的另一层问题。

参考资料

Tobin A. Driscoll、Richard J. Braun，Fundamentals of Numerical Computation，§6.2 Euler’s method，欧拉推进、截断误差及误差传播。
同书 §6.4 Runge–Kutta methods，显式中点法与经典 RK4 的级数。
同作者教材 Absolute stability，线性试验方程、前向和后向欧拉稳定区域。