数值解

Euler 与 RK4：步长、误差与 Taylor 展开的关系。

约 6 分钟

从一个问题开始

问题

单摆方程 $\theta''+\frac gL\sin\theta=0$ 没有初等闭式解。可我们仍然想知道 $\theta(t)$ 是什么样，想画出它的曲线。

我们需要一种只用加减乘除就能算出近似解的方法，而且必须能估计误差、知道什么时候它会失效。

观察

选"单摆"： $\theta''+\gamma\theta'+\frac gL\sin\theta=0$ 。这个方程没有初等闭式解。但实验照样画出了曲线：它是一步一步算出来的。把步长 $h$ 调大，红色的 Euler 曲线开始偏离蓝色的 RK4，甚至发散。回到有精确解的"指数增长"，看右下角两种方法的终点误差随 $h$ 如何变化： $h$ 减半，Euler 误差约减半，RK4 误差约减到 $1/16$ 。

交互实验常微分方程数值解

方程

\theta'' + \gamma\theta' + \tfrac{g}{L}\sin\theta = 0

时间序列

相图

g/L2.00

γ0.10

初值 x(0)2.40

初值 x′(0)0.00

步长 h0.250 · 60 步数

方法 Euler RK4此参数下没有可用的闭式解，只显示数值解。

猜想

$y'=f(t,y)$ 告诉我们当前的斜率。从 $(t_n,y_n)$ 出发，沿这个斜率走一小步 $h$ ，得到下一个近似点。步长越小越准，但误差如何随 $h$ 缩小，取决于我们对"走一步"的近似有多精细：这应该能用 Taylor 展开算清楚。

定义

定义 4.1Euler 法

$y_{n+1}=y_n+h\,f(t_n,y_n),\qquad t_{n+1}=t_n+h.$

定义 4.2局部截断误差与全局误差

局部截断误差是从精确值 $y(t_n)$ 出发走一步后与 $y(t_{n+1})$ 的差。全局误差是 $|y_N-y(T)|$ ，其中 $N=T/h$ 步。若全局误差为 $O(h^p)$ ，称方法是 $p$ 阶的。

推导：Euler 法的阶

命题 4.1Euler 法是一阶方法

局部截断误差为 $O(h^2)$ ，全局误差为 $O(h)$ 。

证明

Taylor 展开精确解： $y(t_n+h)=y(t_n)+hy'(t_n)+\tfrac{h^2}2y''(\xi)=y(t_n)+hf(t_n,y(t_n))+O(h^2)$ 。Euler 一步正好是前两项，所以局部误差 $O(h^2)$ 。

全局误差：设 $e_n=y_n-y(t_n)$ ， $f$ 关于 $y$ 是 $L$ -Lipschitz 的，局部误差不超过 $Ch^2$ 。则 $|e_{n+1}|\le|e_n|+h|f(t_n,y_n)-f(t_n,y(t_n))|+Ch^2\le(1+hL)|e_n|+Ch^2.$ 迭代得 $|e_N|\le Ch^2\sum_{k=0}^{N-1}(1+hL)^k\le Ch^2\cdot\dfrac{(1+hL)^N-1}{hL}\le\dfrac{Ch}{L}\bigl(e^{LT}-1\bigr)$ ，用了 $(1+hL)^N\le e^{hLN}=e^{LT}$ 。所以全局误差 $O(h)$ ： $N$ 步各犯 $O(h^2)$ 的错，累积成 $O(h)$ 。

更高阶：Runge–Kutta

Euler 法只用了区间左端的斜率。更好的做法是在区间内多采几个斜率，加权平均。

定义 4.3经典四阶 Runge–Kutta（RK4）

\begin{aligned} k_1&=f(t_n,\,y_n),\\ k_2&=f\!\left(t_n+\tfrac h2,\ y_n+\tfrac h2k_1\right),\\ k_3&=f\!\left(t_n+\tfrac h2,\ y_n+\tfrac h2k_2\right),\\ k_4&=f\!\left(t_n+h,\ y_n+hk_3\right),\\ y_{n+1}&=y_n+\frac h6\,(k_1+2k_2+2k_3+k_4). \end{aligned}

定理 4.2RK4 是四阶方法

局部截断误差 $O(h^5)$ ，全局误差 $O(h^4)$ 。

证明

完整证明要把 $k_1,\dots,k_4$ 全部 Taylor 展开到 $h^4$ 并与 $y(t_n+h)$ 的展开逐项比较，篇幅较长。这里证明一个足以说明思想的特殊情形： $f=f(t)$ 不依赖 $y$ 。此时 $k_1=f(t_n)$ ， $k_2=k_3=f(t_n+\tfrac h2)$ ， $k_4=f(t_n+h)$ ， $y_{n+1}-y_n=\frac h6\Bigl[f(t_n)+4f\bigl(t_n+\tfrac h2\bigr)+f(t_n+h)\Bigr],$ 这正是 Simpson 求积公式，而 $y(t_{n+1})-y(t_n)=\int_{t_n}^{t_{n+1}}f$ 。Simpson 公式对三次多项式精确，误差为 $-\frac{h^5}{2880}f^{(4)}(\xi)$ ，即 $O(h^5)$ 。一般情形的结论相同，只是展开更繁琐。全局误差从局部误差按 Euler 法证明中同样的累积论证得到 $O(h^4)$ 。

这就是实验中" $h$ 减半，误差减到 $1/16$ "的来源： $(h/2)^4=h^4/16$ 。

定理：为什么大步长会发散

命题 4.3Euler 法的稳定性

对 $y'=\lambda y$ （ $\lambda<0$ ），Euler 法给出 $y_n=(1+h\lambda)^ny_0$ 。真解衰减到 $0$ ，而数值解衰减当且仅当 $|1+h\lambda|<1$ ，即 $h<\dfrac2{|\lambda|}$ 。

证明

直接代入： $y_{n+1}=y_n+h\lambda y_n=(1+h\lambda)y_n$ 。

步长超过 $2/|\lambda|$ 时 $(1+h\lambda)^n$ 振荡放大，与真解毫无关系。这解释了实验里大步长时 Euler 曲线的发散。这种"刚性"问题需要隐式方法，那是数值分析的内容。

常见错误

步长减半，Euler 法的误差只减半；把 $h$ 调得很小并不能弥补方法的低阶，反而增加计算量与舍入误差。而 $h$ 太大时 Euler 法会发散。精度（阶）与稳定性（步长上限）是两件事，都要考虑。

应用

应用没有闭式解的方程

单摆、三体问题、Lorenz 系统、几乎所有真实的物理模型，都只能数值求解。实验里单摆的曲线是 RK4 以 $h$ 为步长、每步四次求值算出来的。步数 $N=T/h$ ；精度和计算量之间的权衡，是整个科学计算的日常。

注

数值方法不是"近似的数学"。它有自己的定理：阶、收敛、稳定性都是可以精确证明的命题。Euler 法的全局误差界 $\dfrac{Ch}L(e^{LT}-1)$ 里的每个字母都有明确含义，而且它诚实地告诉你误差随 $T$ 指数增长：预报越远越不准。

练习

对

y'=y,\ y(0)=1

用步长

h=0.5

的 Euler 法走两步，得到的

y(1)

近似值是多少？

RK4 的全局误差是

O(h^4)

。把步长减半，误差大约变为原来的几分之一？