数值方法：沿着变化率，一步步算出解

上一章里，我们把咖啡的降温、水槽的流入流出、种群的增长都写成了方程。模型一旦写好，接下来通常会问一个很具体的问题：十分钟后温度是多少？某个时刻水槽里有多少盐？前面学过的分离变量和积分因子能够回答一部分问题，可如果环境随时间变化得很复杂，或者变化率来自一张测量表，我们可能根本写不出方便计算的解公式。

这时先别把“没有求出公式”理解成“什么也算不了”。假设你已经知道当前温度，也知道它此刻每分钟下降多少度，那么只向前看很短的一段时间，下一刻的温度还是可以估出来的。到了下一刻，重新查看变化率，再走一小步。数值方法就是把这个朴素动作组织成一套可以重复的计算规则。

我们会从只看一次斜率的 Euler 方法出发，逐步改进对“一步之内究竟变化了多少”的估计。过程中始终保留一个问题：算出来的小数看上去很精细，凭什么相信它？误差、步长和稳定性，都是为了回答这个问题。

先把要计算的东西说清楚

想象一小块金属放进逐渐升温的环境中。金属起初比环境热，会先冷却；环境继续升温后，它又可能转而升温。我们用一个教学模型追踪这个过程：金属内部温度近似均匀，环境温度在观察的两分钟内按恒定速率升高，传热系数保持不变。

设时间 $t$ 以分钟计，金属温度为 $T(t)$，初始温度为 $21\,{}^\circ\mathrm C$，环境温度为 $T_a(t)=20+t$。这里省略书写的升温速率是 $1\,{}^\circ\mathrm C/\mathrm{min}$。冷却常数取 $k=1\,\mathrm{min}^{-1}$，于是“温度变化率与物体和环境的温差相反”给出

为让计算简洁，把相对 $20\,{}^\circ\mathrm C$ 的温度数值记成 $y=T-20$。在上述单位下，方程变为

这条方程允许我们用积分因子求精确解。两边乘 $e^t$，就有 $(e^ty)'=te^t$。因为 $$，所以

代回去检查也很直接：$y'=1-2e^{-t}$，而 $t-y=1-2e^{-t}$；在 ，解确实等于 。这次特意选一个能算出精确解的模型，是为了让数值方法有一个可以核对的答案。等计算规则弄明白了，同样的规则还能用于没有方便解公式的问题。

现在把时间轴切成小段：

$h$ 是每一步跨过的时间，叫作步长；$y_n$ 是我们算出的数，期望它接近精确值 $y(t_n)$。下标 $n$ 数的是“走了几步”，不是求导次数。本章固定 $t_{}$，若用 算到 ，就需要四步。

一般初值问题写成 $y'=f(t,y)$。算法已知的是当前位置 $(t_n,y_n)$ 和右端函数 ，并不知道下一点的真实高度。方向场在这个位置给出斜率 ，数值方法就从这里获取前进的信息。

还有一个容易漏掉的前提：我们计算的区间必须落在解存在的范围内，各个采样点也要落在 $f$ 有定义的区域里。若右端出现 $1/y$，算到 $y=0$ 附近就不能照常跨过去；若真实解在有限时间发散，一串有限的小数也不能证明解已经穿过了那个时刻。

Euler 方法：先把短短一段当成切线

如果此刻变化率是每分钟 $-1$ 度，那么半分钟内大约下降 $0.5$ 度。这个“大约”藏着 Euler 方法的全部想法：在这一小步里，暂时假设变化率不变。

从点 $(t_n,y_n)$ 出发，斜率为 $f(t_n,y_n)$ 的直线是

令 $t=t_n+h$，就得到下一点：

先算斜率，再乘时间，最后加回原高度。 这三个量不能混淆。若 $y$ 是温度，$f$ 的单位是温度每分钟，只有 $hf$ 才能与温度相加。

严格说来，第一步的直线是原初值问题解曲线的切线；后面因为 $y_n$ 已经有偏差，它是经过数值点的那条局部解曲线的切线，不一定还贴在原来的精确解上。Euler 方法每一步都认真遵守当前位置的方向场，但“当前位置”本身可能已经偏了。这也是误差会继续传播的原因。

把四步真正算完

回到 $y'=t-y$，$y(0)=1$，取 $h=0.5$。起点斜率为 $$，所以

来到数值点 $(0.5,0.5)$，斜率变成 $0.5-0.5=0$，于是

这里算出的“这一段没有变化”并不等于金属真的整整半分钟保持恒温。它只是说，我们从这个近似起点取到的斜率恰好是零，并把它沿用了一整步。再到 $(1,0.5)$，斜率为 $0.5$；到 $(1.5,0.75)$，斜率为 $0.75$，所以最后两步为

把每一步放在一起，计算过程就很容易复核：

Euler 给出的两分钟后温度是 $20+y_4=21.125\,{}^\circ\mathrm C$。精确解则给出

误差约为 $0.145671\,{}^\circ\mathrm C$。这一误差在某个应用中能不能接受，要由应用所需的精度决定；这里先关注它为什么出现，以及怎样减小。

局部误差和全局误差，问的是两件事

先做一个思想实验。假设某一步开始前，有人把你已经算偏的数改回真实值 $y(t_n)$，然后只让你用 Euler 走一步。这一步留下的偏差叫作局部一步截断误差：

如果解足够光滑，Taylor 展开告诉我们

由于 $y'=f(t,y)$，Euler 保留的恰好是前两项，剩下的从 $h^2$ 开始。因此局部一步误差为 $O(h^2)$。这里的 表示在足够小的步长范围内，误差绝对值可以被一个与 无关的常数乘 控制；它不等于“误差恰好是 ”。曲线弯得越厉害， 相关的系数也可能越大。

真正运行算法时，可没有人每步帮我们恢复真实值。记

这叫网格点上的全局误差。把真实的一步和数值的一步相减，会得到

这个式子值得读成一句话：下一步的误差，包含旧误差、旧误差引起的斜率偏差，以及这一步新添的截断误差。全局误差不是把“最后一步的局部误差”换个名字。

图：局部误差只看一步，全局误差会带着前面每一步的偏差继续传播。

为什么 Euler 的全局误差只有一阶

设在所研究区域中，$f$ 对 $y$ 满足 Lipschitz 条件：

这和前面唯一性讨论中的条件相连：高度发生一点偏移，斜率不能随之毫无限制地变化。如果还知道 $|\tau_{n+1}|\le Ch^2$，就有

从精确初值 $e_0=0$ 开始，把这个不等式逐步展开。在固定总时长 $b-t_0=Nh$ 内，若 $L>0$，可得到

若 $L=0$，直接相加得到 $|e_N|\le C(b-t_0)h$。因此 Euler 的全局误差是 。直觉上，每步产生约 的误差，但走完固定区间需要约 步；在误差传播可控制的前提下，最终就剩下一个 的量级。

所以，步长减半后，单独一步的误差通常约为原来的四分之一，而同一终点的全局误差通常约为原来的一半。两句话都对，因为比较的对象不同。有些场合会把局部一步误差再除以 $h$ 后命名为“局部截断误差”；本章始终采用未除以 $h$ 的定义。

这些结论讨论的是固定时间区间、足够光滑的方程和足够小的步长。不能把区间同时拉长，再要求误差一定按原来的比例下降。参数特别大、解靠近奇点或右端不光滑，也都会改变实际观察到的效果。

下方实验使用面板中列出的初值问题。先核对方程、初值和终点，再比较步长；若它们与本文温度算例不同，就应单独比较该问题的误差，不能直接套用本文表格中的数值。

Heun 方法：先预测终点，再修正整步变化

Euler 的主要问题是，一整步都只听起点斜率。如果我们能兼顾终点附近的变化率，对平均变化的估计就可能更好。

从积分的角度，一步真实变化满足

Euler 用“区间长度乘左端高度”近似这个积分。另一个自然选择是用两端斜率的平均值乘 $h$，有点像梯形面积。但真正的终点高度还不知道，终点斜率也就算不出来。解决办法是：先借用一次 Euler，做出临时终点；再用这个临时终点提供第二次斜率信息。

把过程写开：

本课程把这套端点平均的预测校正公式称为 Heun 方法，也称改进 Euler 或显式梯形法。数值方法的名称有时存在不同约定，认清上面这三个动作，比只记一个名字可靠。

临时值 $\widetilde y_{n+1}$ 只负责帮我们取斜率，校正后的 $y_{n+1}$ 才交给下一步。最后的增量要加在旧值 上，不能加在预测值上，否则会把前进的距离重复计算。这里也只做一次校正；它不是要反复迭代求解一个隐式方程。

用同一个模型算两步

仍取 $y'=t-y$、$y_0=1$、$h=0.5$。第一步的起点斜率是 ，Euler 预测 。在预测点 ，终点斜率为 。两者平均是 ，所以

这比 Euler 的 $0.5$ 更靠近精确值 $y(0.5)\approx0.713061$。第二步必须从校正后的 $(0.5,0.75)$ 开始：

因此

继续同样的步骤，可以得到完整记录：

为什么多取一次斜率就能提高阶数？把第二次斜率在起点附近展开，会有

这里右边的 $f$ 和偏导数都在起点计算，而沿解曲线 $y''=f_t+f_yf$。因此 Heun 的一步更新恰好包含

也就是精确解 Taylor 展开的前三项。局部一步误差从 $O(h^3)$ 开始，在适当的光滑性和误差传播条件下，全局误差是 $O(h^2)$。它能改善 Euler，是因为确实补上了之前遗漏的一层变化，不只是因为“多算了一遍”。

中点法也算两次斜率，但不是同一个方法

既然目的是估计平均斜率，除了平均两端，也可以去区间中部看一眼。显式中点法先沿起点斜率走半步，估计中点高度，再用中点斜率完成整步：

Heun 用两端斜率的平均，中点法最终只用估计中点处的斜率。两者每步都求两次 $f$，也都是二阶方法，但公式不能混用。

对刚才的线性算例，中点法第一步得到 $k_1=-1$，估计中点为 $(0.25,0.75)$，所以 $k_m=-0.5$，最终 。它恰好和 Heun 相同；实际上这个算例的每一步都会相同。千万别因此认定两个方法本来就一样。

换成一个非线性的短算例就能看清差别。某个无量纲增长模型中，当前量越大，增长率按它的平方增加，写成

这里只把它用于很短的时间，因为其精确解 $y=1/(1-t)$ 在 $t=1$ 发散。取 $h=0.2$：Heun 预测终点为 $1.2$，再取斜率 ${1.2}^2=$，于是

中点法预测中点为 $1.1$，中点斜率为 $1.21$，所以

精确值为 $y(0.2)=1.25$。两种结果都比 Euler 的 $1.2$ 接近，但彼此不同。“同为二阶”只约束步长趋小时的误差量级，并不保证每一个步长、每一道方程都得到一样的数。

RK4：四次斜率采样怎样配合

如果想进一步提高精度，又不想反复对复杂的 $f$ 求高阶导数，就可以在一步内部安排更多次斜率计算。经典四阶 Runge–Kutta 方法通常简称 RK4。它把四次计算组织成以下顺序：

第一次就在起点；第二次去用 $k_1$ 预测的中点；第三次仍在中点时间，却用 $k_2$ 重新预测高度；第四次去用 $k_3$ 预测的终点。四个采样点不必落在精确曲线上，也不是连续向前走的四小步。每次构造高度，都从同一个旧值 $$ 重新出发。

两个中点斜率虽然对应相同时间，高度通常不同，所以不能省掉 $k_3$，也不能把它直接写成 $k_2$。我们这里的 $k_i$ 都是斜率，更新时要乘步长；有的记法把步长提前吸收到 $$ 中，阅读其他公式时要先辨清这一点。

权重 $1:2:2:1$ 与这些采样位置共同配合，使一步展开与精确解匹配到 $h^4$ 项。这是一个需要验证的代数性质，不是“中点天然比端点重要”就能推出的规则；随便换几个权重，一般就不再是四阶。在足够光滑的条件下，RK4 的局部一步误差为 $O(h^5)$，全局误差为 $$。

第一小步，四次斜率全部展开

对温度算例，$y_0=1$、$t_0=0$、$h=0.5$。四个斜率依次为

所以

精确值约为 $0.7130613194$，第一步的绝对误差约为 $0.0004803472$。这个结果确实明显改善了，不过只看第一步还不够。下一步要用这个近似值重算一整组 $k_1,k_2,k_3,k_4$，前一步的斜率不能原封不动沿用。

用一个容易验算的模型看“四阶”的意思

设一个无量纲增长过程满足 $u'=u$。若某一步起点值为 $u_n$，把它代入 RK4 的四次计算，整理后得到

精确地前进同样的时间，则应乘 $e^h$。指数展开的前五项刚好就是括号里的多项式，差异从 $h^5$ 项开始。这给出了四阶性质的一个可直接核对的例子。一般非线性方程还要核对更多项，不能只凭这个线性例子就当成完整证明。

比较方法时，要同时看精度和计算量

三种方法现在都已经算得出来。对同一个温度模型，以同样的 $h=0.5$ 走四步，结果如下。表格只在展示时四舍五入，后续计算保留更高精度。

Euler 低估，Heun 和 RK4 在这个算例中高估。不要把这当成方法的一般属性，误差方向取决于方程和所走的区间。这里三个方法都抓住了“先降后升”的大致趋势，但转折位置和数值精度不同。精确解的转折满足 $y'=1-2e^{-t}=0$，在 $t=\ln2\approx0.6931$；粗网格只记录少数时刻，不能准确确定最低点出现的时间。

再把步长减半，始终比较同一个终点 $t=2$ 的绝对误差：

Euler 的误差比逐渐接近 $2$，Heun 接近 $4$，RK4 接近 $16$，对应全局阶数 $1、2、4$。第一组比值未必就等于这些数，因为步长还没有足够小，高阶余项仍可能明显影响误差。

同步长比较时，RK4 花的计算量更多。如果总预算都是 $16$ 次函数计算，Euler 可以在 $[0,2]$ 上走 $16$ 步，Heun 走 $8$ 步，RK4 走 $4$ 步。根据上表，三者误差分别约为 $0.03453639、0.00688519、0.00042897$，这个算例中 RK4 在相同预算下仍然占优。结论来自实际比较，不是由“四阶”这个名字自动保证每种情况都最好。

如果没有精确解，怎样估计还差多少

实际问题中，表格最后一列往往没有办法填写。一个常用办法是用步长 $h$ 和 $h/2$ 各算一次，在相同的时刻比较两次结果。

假设某个 $p$ 阶方法已经进入稳定的渐近范围，终点近似满足

于是细网格结果的主导误差可以用

来估计。注意分母是 $2^p-1$，这个估计针对的是较细步长结果。对 Euler、Heun、RK4，分母分别是 $1、3、15$。

还可以再算一次 $h/4$，观察相邻差值的比是否接近 $2^p$。若差值根本不按预期减小，就该检查是否尚未进入小步长范围、跨过了不光滑点，或发生了稳定性问题。两次结果碰巧接近不等于已经获得严格误差保证，多显示几位小数同样不等于多知道几位准确数字。

对复杂轨迹，可以在变化剧烈处取小步长、变化缓慢处取大步长，这就是自适应步长的基本方向。实际算法常用一对不同阶数的近似来估计局部误差，再与绝对和相对容差比较。本章不实现完整控制器，但要记住：光看当前斜率大不大不足以决定步长，斜率变化、误差传播和稳定性都可能提出额外限制。

步长太大，冷却竟然会算成爆炸

误差大有时只是最后多差一点点；更严重的情况是，算法连变化趋势都算反了。

把物体放在恒温环境中，用 $y$ 表示物体与环境的温差。若热量交换遵循此前的假设，温差满足

精确解 $y(t)=y_0e^{-at}$ 始终趋近零。若最初温差为正，它会从正的一侧单调下降，不会越过零，更不会忽大忽小。

Euler 更新却是

从而

数值解是否衰减，全看每一步乘上的因子 $1-ah$。要让其绝对值逐步缩小，需要

这个条件叫作该测试方程上 Euler 的绝对稳定条件。$a$ 的单位是时间的倒数，因此 $ah$ 没有单位，条件在量纲上也合理。

如果 $0<ah<1$，数值温差保持原来的符号并衰减；如果 $ah=1$，算法一步变成零，随后停在那里，但真实温差在有限时刻还没有精确变成零。如果 $1<ah<2$，数值解正负交替、幅度逐渐变小，虽然最终趋零，却已经出现真实模型没有的振荡。当 $ah=2$，因子为 ，振幅不再减小；当 ，绝对值反而越乘越大。

例如取 $a=4\,\mathrm{min}^{-1}$、$h=0.75\,\mathrm{min}$、$y_0=10\,{}^\circ\mathrm C$，Euler 的因子是 ，算出的温差依次为

精确温差在第一步之后却只有 $10e^{-3}\approx0.497871$ 度。模型在迅速冷却，算法在制造越来越剧烈的冷热交替。这种错误即使用精确分数计算也会发生，不能归咎于计算机的小数舍入。

高阶显式方法也有稳定范围。Heun 用于同一测试方程时，每步乘上 $1-ah+(ah)^2/2$；它虽然是二阶方法，衰减所需的区间仍为 $0<ah<2$。RK4 的因子则为

仍须检查 $|R(-ah)|<1$。例如 $ah=3$ 时，因子为 $1.375$，RK4 同样把衰减算成增长。提高阶数不能代替稳定性检查。

刚性：你关心的过程很慢，算法却被迫走得很小

有些系统同时存在很快的调整和很慢的变化。金属迅速追随一个缓慢变化的环境，就能提供这种情形的入口。为了突出时间尺度，考虑已经无量纲化的模型

直接代入可知精确解是 $y=\cos t$，它变化得并不急。可是令偏离量 $z=y-\cos t$，方程就给出

真实系统会迅速压下这种偏离；Euler 对两条数值轨迹间的偏差却每步乘以 $1-1000h$。为不放大偏差，它必须取 $h<0.002$。即使我们只想观察缓慢的余弦变化，也不能随意跨大步。

这种“精度看起来允许较大步长，显式方法却因快速衰减成分的稳定性被迫使用很小步长”的情形，是理解刚性的一个实用入口。实际判断还与观察区间、所需精度和所选方法有关，不是看到一个大系数就机械贴标签。

一种后续会用到的思路，是让新时刻的变化率参与更新。例如后退 Euler 对 $y'=-ay$ 写成

对任意 $h>0$，这个衰减因子都在 $0$ 与 $1$ 之间。代价是一般非线性问题里，需要求解关于新值的方程。它提示我们：方法的选择不只是在精度阶数之间做比较，有时还要改变更新的方式。

已经会的解析方法，也能帮数值计算减轻负担

数值方法并不要求把前面学过的积分因子全部放下。若方程里有一部分结构能够精确处理，我们可以先处理它，再近似剩下的部分。

例如一个量在以固定比例衰减，同时受到随时间改变的补充，方程写成

乘积分因子 $e^{at}$，并令 $u=e^{at}y$，就得到

现在 $u$ 的变化率只与时间有关。我们可以用矩形、梯形或更高阶的求积近似来计算 $u$ 的增量，再通过 $y=e^{-at}u$ 换回原量。这样，已知的指数衰减被精确保留下来，需要近似的只是外部补充在时间中的累计作用。

还可以直接在一小段上积分，得到完全准确的关系

这条式子很好理解：原来已有的量经历了 $h$ 时间的衰减；时刻 $s$ 新补入的量只经历从 $s$ 到终点的衰减，所以带着不同的指数权重。它把建模意义和积分因子联系在了一起。

若一步内暂时把补充率视为常数 $g(t_n)$，积分可以直接算完，得到一种更新规则

例如 $a=4$、$g\equiv8$、$y_0=3$，精确解是 $y=2+e^{-4t}$。即使用 ，上式仍给出 ，恰好等于精确值，因为恒定补充的积分没有再作近似。普通 Euler 在这个步长下却给出 ，并在后面围绕平衡值 振荡发散。

这个例子不是说任何变量替换都必定改善结果。变换可能使新的右端变化得更快，也可能使中间计算出现很大的指数。我们仍须比较步长和误差。它说明的是：解析分析和数值计算可以配合；先认出方程中能够精确保留的结构，有时比直接缩小步长更有效。

把数值结果放回方向场和模型里

经过前面的计算，我们已经有了检查结果的几个具体抓手。先看方程告诉了什么：在哪些区域上升，在哪些区域下降，平衡值在哪里；再看数值轨迹有没有遵守这些结构。

温度算例中，$y'=t-y$，因此 $y>t$ 时金属比环境热，温度下降；$y<t$ 时金属比环境冷，温度上升；直线 $y=t$ 上的瞬时斜率为零。精确解会在 碰到这条零斜率线，然后从冷却转为升温。方向场使“先降后升”有了方程上的理由，也提醒我们粗步长可能把转折位置算偏。

如果模型要求浓度非负，数值结果却突然出现负浓度，先检查右端、单位和步长，不要马上为负值编一个现实解释。若外部输入在某个时刻突然切换，应把计算区间在切换点分开，让算法在每一侧使用正确的方程；跨越不光滑点的一整步，通常不能继续套用光滑情形的高阶误差结论。

实践中可以先用中等步长观察趋势，再缩小步长检查同一组输出时刻。如果能写出精确解，就直接比较；如果不能，就结合步长加密、不同方法以及模型约束来检查。输出表还应说明时间单位、使用的方法和步长，否则别人很难判断那串数字是怎样得到的。

最后区分三种来源不同的偏差。截断误差来自用有限次采样近似连续变化；舍入误差来自有限位数运算；模型误差来自传热系数恒定、混合充分等假设与真实对象的差别。把步长缩小主要针对第一种。即使方程被解得极准，也不会自动修正一个选错的模型。

练习：每个近似值都要有来路

起点斜率怎样变成下一点

一个无量纲初值问题为 $y'=t+y$、$y(0)=1$。用 Euler 方法取 $h=0.25$ 算到 $$，并用精确解检查误差。

同一个起点，两种二阶方法

对 $y'=y^2$、$y(0)=1$，分别用 Heun 和显式中点法，以 $h=0.1$ 做一步。两者为什么不同？

RK4 的四次计算不能省略

对无量纲衰减模型 $y'=-2y$、$y(0)=1$，用 RK4 取 $h=0.5$ 做一步，列出四个斜率并比较精确值。

四分之一和二分之一，哪个说法对

有人说 Euler 的局部一步误差是 $O(h^2)$，因此把步长减半后，固定终点的误差也应减到四分之一。解释问题出在哪里。若某二阶方法在步长 $h$、$h/2$ 时的终点值为 $2.408$ 和 $2.402$，在渐近假设下估计细网格结果的误差。

冷却模型什么时候会被算成振荡

温差满足 $y'=-5y$，$t$ 以分钟计，初始温差为正。用 Euler 时，$h=0.1、0.3、0.4、0.5$ 分钟分别会出现什么？

一串有限数字能不能越过真实解的终点

对 $y'=y^2$、$y(0)=1$，某程序使用固定步长，输出了 $t=1.2$ 处的一个有限值，并把它标成“原初值问题的近似解”。这个标注合理吗？应该怎样检查？

经历这一章，我们已经能在没有显式公式时，把变化规律变成一串可检查的近似值。本章主要计算的是单个一阶方程，只需记录一个状态量。接下来面对弹簧上的物体，同一个位置可能正在向左运动，也可能正在向右运动；只给位置，未来并没有确定下来。下一章会从这个问题进入二阶线性方程，说明为什么除了初始位置，还需要初始速度，以及两个初始条件怎样选出一条运动轨迹。

用积分因子解 $y'-y=t$，得到 $y=Ce^t-t-1$。初值给出 $C=2$，所以 $y(0.5)=2e^{0.5}-1.5\approx1.797443$，Euler 的绝对误差约为 $0.172443$。这里差的是高度，不是斜率。

中点法预测中点为 $1.05$，取中点斜率 $1.1025$，于是

精确值为 $1/0.9\approx1.111111$。两种方法采样的位置和最后组合斜率的方式不同；同属二阶，并不意味着误差系数相同。

所以

精确值为 $e^{-1}\approx0.367879$，绝对误差约 $0.007121$。第二、第三次的时间相同，但预测高度分别为 $0.5$ 和 $0.75$，斜率当然不能直接视为相同。