同一个衰减方程，用一种步长计算，结果越来越小；把步长加倍，结果却开始正负交替，幅度越跳越大。方程没改，初值也没改。问题出在每次向前走时，算法怎样处理已经积累的偏差。

这一章把“走一步”拆开来看。你会亲手算出 Euler、中点法和 RK4 的一步，分清局部误差与全局误差，再用稳定域判断哪些步长会把衰减算成增长。遇到刚性问题时，我们还要衡量：多花力气解一个隐式方程，能否省下大量过小的显式步。

从真实解出发走一步

考虑初值问题

$y(t)$ 表示真实解，$Y_n$ 表示程序在 $t_n=nh$ 算出的值。两者不要混用。对向量方程，下面的绝对值换成所选范数即可；需要的光滑性和 Lipschitz 条件仍须在有关区域成立。

显式 Euler 从当前斜率作预测：

更一般地，把一个单步方法记为

想看这一步本身有多准，可以暂时把起点换成真实值，定义一步缺陷

这个量不含前面各步的误差。对 Euler，Taylor 展开给出

所以 $d_{n+1}=h^2y''(t_n)/2+O(h^3)$。有些书把 称为局部截断误差，那时 Euler 的局部截断误差写成 。本章始终用未除以 的一步缺陷，因而是 ；看到阶数相差一，不妨先核对定义。

真正关心的是 $e_n=y(t_n)-Y_n$。它包含之前的偏差，以及这些偏差后来受到的放大。假定在固定区间 $[0,T]$ 内，真实解与数值路径都留在一个区域中，并且对充分小的 ，有与 无关的常数 ，使

Euler 的第二个条件可由 $f$ 对 $y$ 的 Lipschitz 界推出；一般 RK 方法要检查自己的整个更新映射。现在把误差拆成两项：

由此得到

连续代入这个不等式，每一份早先的误差都带上若干次放大因子。当 $L>0$ 时，几何级数给出

若 $L=0$，直接求和得到 $|e_n|\le|e_0|+Ct_nh^p$。因此在初值精确、条件一致成立的情况下，一步缺陷 导出固定有限时间内的全局误差 。

这里没有承诺误差永远不增长，也没有让 $T$ 随意趋于无穷。把“算法的阶数”“有限时间的误差传播”“固定步长的长期稳定性”分开，后面就不会互相矛盾。

把真实终点与数值终点之间的差拆成两段：从真实起点计算也会留下本步缺陷，原来起点中的误差则通过更新映射继续传播。框的高低只表示拆分关系，不代表数值大小。

用中途的斜率修正方向

Euler 一整步都采用左端斜率。显式中点法多问一次：沿着这条斜率走半步后，斜率会是什么？

本章的 $k_i$ 都是斜率，没有预先乘 $h$。中间预测值 $Y_n+hk_1/2$ 也不是真实的半步解，只是为了取得更合适的斜率。

为什么这样能提高一阶？从真实起点作多元 Taylor 展开，所有导数都在 $(t_n,y(t_n))$ 计算：

而沿着真实解求导，$y''=f_t+f_yf$。于是更新值包含

恰好匹配真实解到二次项；一步缺陷是 $O(h^3)$，在上一节的传播条件下，全局是二阶。

这个推理还能帮我们设计另一种二阶方法。设

展开后，一次项的系数是 $b_1+b_2$，二次项的系数是 $b_2\alpha$。只要

就能匹配到二阶。例如 $\alpha=1$ 时，两个斜率各取一半；它与显式中点法的节点不同，却同为二阶。以后比较算法，不能只凭“计算了几次斜率”判断精度。

RK4 的四次计算各有位置

经典四阶 Runge–Kutta 方法的完整规则是

第二、三次虽然都在半步时刻，却通常用不同的 $Y$ 值。第四次也不是在已经算好的 $Y_{n+1}$ 上取斜率，因为最终值此时还没求出。

把一整步算出来。 取 $y'=t+y$、$y(0)=1$、$h=0.2$。真实解为 $$，代回方程和初值即可核对。RK4 的四次取样分别为

所以

真实值约为 $1.2428055163$。同一步的 Euler 值为 $1.2$，显式中点值为 $1.24$。这次比较说明了具体算例的误差大小；单独一个例子还不能证明 RK4 对一般光滑方程都是四阶。

中间两次取样在同一个时刻，却使用不同的状态。箭头表示前一斜率参与下一取样状态的计算，并不表示两者数值相等；最终更新值还要按权重合成。

四阶来自哪些匹配

下面给出一般光滑方程的核对过程。符号会稍密一些，但每一项都有出处。把时间也作为未知量，写成 $X=(t,y)$、$X'=F(X)=(1,f(t,y))$，就能统一处理显含时间的方程。

在当前点记 $J=DF$ 为一阶导数，$H=D^2F$ 为双线性二阶导数，$\mathcal T=D^3F$ 为三线性三阶导数。例如 是矩阵作用于向量， 是把 放进二阶导数的两个位置； 表示 ，并不是对 再求一次导数。

沿真实解使用链式法则：

例如，对 $J^2F$ 求导，会得到 $H(F,JF)+JH(F,F)+J^3F$；对 求导，则得到 。这解释了系数 3 从哪里来。

设 RK 的阶段为 $K_i=F(X+h\sum_j a_{ij}K_j)$，最终值为 ；这里大写 是向量阶段，与三阶导数 不同。令 、。对每个阶段的增量作 Taylor 展开，逐次代入前面已知的阶段，可得

其中 $c^2,c^3$ 表示逐分量平方、立方。展开中的交叉项 $H(F,JF)$，来自二阶 Taylor 项里一份 $hF$ 与一份 $h^2JF$ 的乘积。

把上式乘上 $hb_i$ 并求和，再与真实解的 Taylor 式比较，需要以下八个系数相符：

RK4 的系数为

直接乘得

前三个单纯幂次的加权和分别为 $b^Tc=1/2$、$b^Tc^2=1/3$、$b^Tc^3=1/4$；其余四项为

加上 $\sum b_i=1$，全部匹配。若有关导数足够光滑且在所用邻域有一致界，剩余的一步缺陷是 $O(h^5)$；再结合单步映射的 Lipschitz 界，得到有限时间内的全局四阶。只检验 $y'=\lambda y$ 时， 都为零，会漏掉上面若干条件，因此线性测试不能代替这份一般论证。

在同一个终点比较

仍算 $y'=t+y$ 到 $T=1$，把步数从 5 加到 10、20。以下是终点绝对误差，末位作了舍入：

步长减半后，误差比分别向 $2,4,16$ 靠近，不要求较粗网格已经恰好等于这些数。计算成本也要记录：这三种显式方法每步分别调用右端函数 1、2、4 次。相同步数不等于相同成本。

在实验里逐步显示四个取样点，留意 $k_2$ 与 $k_3$ 的位置差别；随后切换方法，比较完整轨迹和终点误差。图上连接数值节点的线段用于显示路径，不额外提供 RK4 阶的连续插值。

没有真实解时，三张网格能告诉你什么

设同一终点的近似值满足 $Y_h=y+Ch^p+o(h^p)$，且 $C\mathrm{\ne }$。减去细步结果，主要误差项成为

因此，细步结果本身的有符号误差约为

若想把细步结果向真实值修正，要加的是相反数：

这两个式子的减法顺序不能混。比如已知二阶计算给出 $Y_h=1.04$、$Y_{h/2}=1.01$，我们估计细步结果偏高 $0.01$，修正值为 $$。

增加第三张网格，可以检查

这只是进入渐近区间的证据。若最高阶系数恰好为零、相邻差受到舍入影响、步长越过稳定范围，或者解本身不够光滑，比值都可能偏离预期。两份近似互相接近，也可能是共同漏掉了同一变化；它们不能自动组成严格误差上界。

三张网格的相邻差比为 4，与二阶主误差项相符。最细结果被估计为偏高 0.005，所以修正时应减去 0.005；这是渐近估计，不是严格上界。

把衰减算成增长：稳定域

用测试方程 $y'=\lambda y$，考察固定步长的长期表现。令 $z=h\lambda$，每一步若写成 $Y_{n+1}=R(z)Y_n$，则 。 叫作放大因子。

当 $|R(z)|\le1$ 时，这个标量序列对任意初值保持有界，我们把这些 $z$ 组成的集合称为绝对稳定域；$|R(z)|<1$ 才会衰减到零。边界上的有界与衰减是两件事。

逐项代入三个显式方法，得到

例如中点法的中途预测为 $(1+z/2)Y_n$，最终值便是 $[1+z(1+z/2)]Y_n$。RK4 同样可从四个阶段展开，不必把放大因子当成另一个要背的公式。

Euler 的稳定域是复平面上 $|z+1|\le1$ 的圆盘，圆心为 $-1$、半径为 1。对于负实数 $\lambda$，若希望数值解也衰减，就须

取 $\lambda=-15$：$h=0.1$ 时 $R=-0.5$，正负交替但幅度缩小；$h=0.2$ 时 $R=-2$，幅度逐步翻倍。精确解在两种情况下都是 。边界 则给出 ，序列始终在 之间循环，虽有界，却没有重现衰减。

Euler 稳定圆盘的圆心为 −1，半径为 1。圆盘内部的 −1.5 对应交替衰减，边界 −2 对应持续振荡，外部 −3 对应交替增长。

如果 $\lambda=i\omega$ 且 $\omega\ne0$，Euler 的 $|1+ih\omega|>1$，任何正步长下的长期幅度都会增长。这不否定它在固定 $T$ 上、 时的一阶收敛：前者固定 让时间趋于无穷，后者固定终点让步长趋于零。

对可对角化的常系数系统，每个特征模态可用这种标量分析。变回原坐标时还要考虑特征向量矩阵的条件数；一般非线性、时变或非正规系统，仅看瞬时特征值并不能得到无条件的整体稳定保证。

刚性使步长受到另一种限制

看一个可以写出真实解的方程：

令 $w=y-\cos t$，便有 $w'=-\kappa w$，所以

当 $\kappa$ 很大，偏离余弦的那部分迅速衰减；之后解变化很慢。但 Euler 传播两条数值路径之间的差，仍要乘 $1-\kappa h$。若 $\kappa=1000$、$h=0.01$，这个因子是 $-9$。即使 $\delta=0$，每步截断产生的偏差也可能被不断放大。

当稳定性要求的步长远小于为了描述所关心变化而需要的步长时，我们称这个计算具有刚性。它与时间区间、精度目标和方法都有关系，不能只凭曲线陡不陡判断。

后向 Euler 改在新时刻取斜率：

对测试方程，这变成 $(1-z)Y_{n+1}=Y_n$，所以

若 $\operatorname{Re}z<0$，则 $|1-z|^2=(1-\operatorname{Re}z)^2+(\operatorname{Im}z)^2>1$，从而 。稳定域包含整个左半平面的方法称为 A 稳定方法，后向 Euler 就是一个例子。沿负实轴让 ，它的放大因子还会趋于零，能强烈抑制快速衰减模态。

对刚才的余弦问题，隐式方程仍然是线性的，整理可得

这个式子可以直接算。不过，稳定只说明相应扰动不会被这一步持续放大；后向 Euler 仍是一阶方法。如果我们关心 $\delta e^{-\kappa t}$ 的短暂过程，大步长可能把它整个跨过去。方法选择要同时回答“会不会失稳”和“有没有分辨想看的变化”。

上方曲线只示意慢变化与快速衰减两种成分，未按同一时间尺度精确作图。下方数值链比较两条离散解之间同一初始扰动的传播，不代表整个余弦解的数值。

实验同时显示稳定域上的 $z$、放大因子和数值轨迹。先在衰减方程里把 Euler 的 $z$ 移到 $-2$ 两侧，再切换余弦刚性问题，看真实解很平缓时数值误差能否仍然增长。复数模式的轨迹显示模长，以免把实部过零误认为解已衰减。

隐式方程也要有停止依据

若 $f$ 非线性，每一步要解

Newton 迭代是

向量形式则解一个以 $I-hJ_f$ 为系数的线性系统，不能逐分量直接相除。第 3、6 章的求根条件和线性求解检查，在这里仍然需要。

例如 $y'=-y^2$、$Y_n=1$、$h=1\mathrm{/}2$，隐式方程为

我们求非负根 $V_*=-1+\sqrt3\approx0.7320508076$。从 $$ 出发，Newton 依次给出 、，继续可接近该根。真实解在这一步末端是 ，所以即便隐式方程解得完全准确，离散误差仍存在。

内层残差 $G(V)$ 和 ODE 的真实误差也不是同一个量。在这个非负区域，$G'(V)=1+V\ge1$，中值定理给出

因此它能控制本步非线性方程的求解误差。一般方程要有相应的逆导数或条件数估计。内层求解留下的误差，会像额外的一步缺陷一样传播；在统一的稳定性估计下，为保持 $p$ 阶全局收敛，可以把每步内层解误差控制到 $O(h^{p+1})$。固定一个与步长无关的粗容差，可能使外层细化很快失去作用。

让程序决定接受还是重算

知道阶数之后，可以让算法试走一步，再决定这一步是否值得保留。这里用显式中点法演示步长加倍比较：从同一个 $(t,Y)$ 出发，一次走 $h$ 得到 $Y_c$；另外连续走两次 $h/2$ 得到 $Y_f$。

从当前初值发出的真实局部解在 $t+h$ 的值记为 $y_*$。光滑情形下，中点法的一步误差可写为

两次半步分别产生约 $C(h/2)^3$ 的误差，前一次误差传播半步的放大为 $1+O(h)$，而第二次的局部系数变化也是 $O(h)$。合起来

于是细结果的误差估计与修正量分别是

这里分母仍为 3，但它来自同一小段里的两次半步；用来预测新步长的误差规模是 $h^3$，不能与固定终点的二阶误差混淆。

设绝对容差为 $\mathrm{atol}>0$，相对容差为 $\mathrm{rtol}\ge0$，定义

若 $\eta\le1$，接受 $Y_f$ 并把时间推进到 $t+h$；若 $\eta>1$，丢弃这次候选，仍从原来的 用较小步长重算。这里接受的是未外推的两半步结果，不能顺手把整个程序叫作四阶方法。

因为估计误差约随 $h^3$ 变化，可用

预测下一次尝试的步长。$0.9$ 留出余量，上下限限制一次调整的幅度；若 $\eta=0$，直接使用允许的最大增长倍数 2。向量问题可逐分量构造尺度，再取各归一化误差的最大值，避免单位较大的分量掩盖其他分量。

一次真实的拒绝。 对 $y'=y$、$Y=1$，尝试 $h=0.2$：中点整步得到 $Y_c=1.22$，两次半步得到 。若 、，则

这一步不接受。新步长约为 $0.1195$，时间仍然是 0，起点仍然是 1。程序不能把已经拒绝的 $1.221025$ 当作下一步的起点。

两个候选都从同一个起点出发。归一化误差估计超过 1 后，本次尝试被丢弃；重算只缩小步长，时间和数值起点均保持原值。

实现还要处理一些很具体的情况：到达终点前截短最后一步；$t+h=t$ 或所需步长低于下限时停止并说明未完成；函数值非有限、尝试次数用尽时也不能假装求解成功。容差需要与数据精度和目标尺度匹配。

非均匀时间取样的示意：某些区段的节点更密。节点密度本身不是误差证明；实际程序依据局部估计、最大步长和失败处理规则决定如何推进。

更要记住，局部误差估计不是全局误差证书。除了之前各步误差会传播，两个候选还可能同时漏掉信息。取 $y'=\sin^2(4\pi t)$、$y(0)=0$，从 0 一步走到 1。整步中点在 $1/2$ 取样，两半步在 取样；三处斜率都为零，所以 ，而真实增量是 。这时估计相等并不意味着真实误差消失。适当限制最大步长、改变起始网格并交叉检查，往往能暴露这样的漏采样，但仍需根据问题结构判断。

在实验里查看每一次尝试的起点、整步与半步候选、$\eta$ 和接受状态。把容差收紧后，不只看最终数字，也看看哪些步被拒绝、计算了多少次右端函数。遇到漏采样预设时，比较修改最大步长前后的结果。

练习：把每一项检查接起来

1. 初值也带着误差。 一个四阶单步方法满足本章的传播条件，但输入初值与真实初值相差 $h$。在固定终点上，能否仍由本章估计保证总误差是 $O(h^4)$？若 $L=0$、一步缺陷上界为 $2h^3$，初值精确，在 $T=$ 的全局误差上界又是多少？

2. 换一组二阶节点。 在两阶段公式中取 $\alpha=2/3$，求二阶所需的 $b_1,b_2$。用它对 $y'=y^2$、 走 ，并与显式中点法的一步比较。

3. 检验一个冒充四阶的方法。 有人取 RK4 原来的 $A,c$，却把四个阶段等权平均，即 $b_i=1/4$。它满足 $\sum b_i=1$、。这就足以保证四阶吗？找出一个不满足的条件。

4. 给修正量定方向。 在同一终点上，三次计算为 $Y_h=2.08$、$Y_{h/2}=2.02$、$Y_{h/4}=2.005$。假定处于单一主误差项支配的渐近区间，估计阶数、最细结果的有符号误差与外推值。这样的数字是否构成严格误差界？

5. 稳定域的边缘。 对 $y'=-20y$、$Y_0=1$，分别用 Euler 的 $h=0.04,0.1,0.12$。写出放大因子并判断长期行为。若改用后向 Euler、，因子是多少？

6. 一个隐式步的内层误差。 对 $y'=-y^2$、$Y_n=1$ 取 $h=1$。写出后向 Euler 方程，求非负根。从 做两轮 Newton。用所得方程残差给出内层解误差上界；再说明它为何不是这一步 ODE 误差的上界。

7. 用成本判断刚性。 对余弦问题取 $\kappa=500$，关心 $0\le t\le10$ 的慢变化。Euler 要衰减快速扰动时，等步长至少需要多少步？如果还要看清非零 $\delta$ 带来的初始快速过程，后向 Euler 的 A 稳定性能否让你任意取大步？

8. 接受不等于误差为零。 自适应中点法从 $y'=y$、$Y=1$ 出发，尝试 $h=0.1$，取 $\mathrm{atol}=10^{-4}$、。算出整步、两半步和 ，判断是否接受。若估计的 超过 1，重算应使用哪个起点？