本文是《微分方程完全入门》系列的第五篇，覆盖ODE数值解法的基本方法：欧拉法、改进欧拉法、四阶龙格-库塔法（RK4），以及自适应步长和刚性方程的概念。阅读本文前需理解一阶ODE的基本形式。本文重点在方法的几何直觉和误差分析，而非纯数值分析视角。

第五篇：数值解法——从欧拉法到 RK4

一、概念引入——为什么要数值解法

在之前四篇中，我们学习了多种解析解法——分离变量、积分因子、特征根法、拉普拉斯变换。

但这些方法都有一个致命的限制：大多数现实中的ODE没有解析解。

考虑一个稍微修改的摆方程——加入非线性项和非正弦驱动力：

\frac{d^2\theta}{dt^2} + \frac{g}{L}\sin\theta = F_0\cos(\omega t)

这个方程——受迫阻尼摆——无法用初等函数表示其解。但它描述的就是真实世界的物理。

数值解法就是回答这个问题的：用计算机逐”小步”逼近解。

二、欧拉法——最直观的起点

2.1 思想

一个一阶 ODE： $\frac{dy}{dt} = f(t, y)$ ，给定初始条件 $y(t_0) = y_0$ 。

在 $(t_0, y_0)$ 处，切线的斜率是 $f(t_0, y_0)$ 。如果我们向前走一小步 $\Delta t$ ：

y_1 = y_0 + f(t_0, y_0) \cdot \Delta t

这就是欧拉法。

2.2 公式

y_{n+1} = y_n + f(t_n, y_n) \cdot h

其中 $h = \Delta t$ 是步长。

2.3 手算例子

用欧拉法求解 $y' = -2y$ ， $y(0) = 1$ ，步长 $h=0.1$ ，计算到 $t=0.5$ 。

$t_n$	$y_n$	$f(t_n,y_n)=-2y_n$	$y_{n+1}=y_n + hf$
0.0	1.0000	-2.0000	0.8000
0.1	0.8000	-1.6000	0.6400
0.2	0.6400	-1.2800	0.5120
0.3	0.5120	-1.0240	0.4096
0.4	0.4096	-0.8192	0.3277
0.5	0.3277

精确解： $y(0.5) = e^{-1} = 0.3679$

误差： $|0.3679 - 0.3277| = 0.0402$ （约 11%）

三、RK4——工程界的标准

3.1 思想

欧拉法只用了区间起点处的斜率，误差 $O(h)$ 。如果能在区间内多取几个点计算斜率，然后加权平均，精度会大幅提升。

RK4 在每一步内计算 4 个斜率：

k_1 = f(t_n, y_n)

k_2 = f(t_n + \frac{h}{2}, y_n + \frac{h}{2}k_1)

k_3 = f(t_n + \frac{h}{2}, y_n + \frac{h}{2}k_2)

k_4 = f(t_n + h, y_n + h k_3)

最终更新：

y_{n+1} = y_n + \frac{h}{6}(k_1 + 2k_2 + 2k_3 + k_4)

3.2 与欧拉法的对比

方法	每步函数计算次数	局部截断误差	步长 $h=0.1$ 时 $t=0.5$ 的误差
欧拉法	1	$O(h^2)$	11%
RK4	4	$O(h^5)$	远小于 0.01%

3.3 RK4的几何直觉

$k_1$ ：起点处切线斜率（欧拉法的值）
$k_2$ ：用 $k_1$ 走到中点后的斜率
$k_3$ ：用 $k_2$ 再算一次中点斜率（修正值）
$k_4$ ：用 $k_3$ 走到终点后的斜率
最后加权平均： $(k_1 + 2k_2 + 2k_3 + k_4)/6$

之所以 $k_2$ 和 $k_3$ 的权重是 2（比 $k_1, k_4$ 大），是因为它们的精度更好——中点的斜率最能代表整个区间的”平均斜率”。

四、数值稳定性和步长选择

4.1 稳定性

如果步长过大，欧拉法可能完全发散。

考虑 $y' = -100y$ ， $y(0) = 1$ 。如果用 $h=0.02$ ：

y_1 = 1 + 0.02 \times (-100) = 1 - 2 = -1

一步就跳到了负值！然后继续发散。

稳定性条件：对于 $y' = \lambda y$ （ $\lambda < 0$ ），欧拉法要求 $h < 2/|\lambda|$ 。

RK4 的稳定性区域更大——最大步长可达约 $2.78/|\lambda|$ 。

4.2 自适应步长

现代 ODE 求解器使用自适应步长——根据局部误差估计自动调节步长。

常用方法：同时计算 RK4 和 RK5，用两者的差值估计误差；如果误差太大，步长减半；误差太小，步长加倍。

MATLAB 的 ode45 使用的就是这种”RK4(5)”方法。

4.3 刚性方程

如果系统中同时存在”很快”和”很慢”的动态：

y' = \begin{bmatrix} -1000 & 0 \\ 0 & -1 \end{bmatrix} y

显式方法必须用小步长满足快速分量的稳定性要求——即使它早已衰减到几乎为零。

刚性方程需要隐式方法（如隐式欧拉法），它的稳定性不受快速分量限制。

在飞行动力学中，六自由度系统的快模态（荷兰滚、短周期）和慢模态（螺旋模态）并存的场景就是典型的刚性方程。

五、在 PX4 仿真中的应用

PX4 SITL 仿真中使用的是 RK4 求解器，步长通常为 1 ms 或 4 ms（与飞控主频匹配）。

仿真类型	步长	方法	原因
PX4 SITL (Gazebo)	1-4 ms	RK4	平衡精度和速度
JSBSim	-	RK4	工程标准
离线模拟	可调	RK4/自适应	精度优先
实时硬件回环	固定小步长	欧拉/RK4	实时约束

六、核心公式速查卡

公式	方法	误差
$y_{n+1} = y_n + h f(t_n, y_n)$	欧拉法	$O(h^2)$
$y_{n+1} = y_n + \frac{h}{2}[f(t_n,y_n) + f(t_{n+1}, y_n+h f_n)]$	改进欧拉	$O(h^3)$
$y_{n+1} = y_n + \frac{h}{6}(k_1+2k_2+2k_3+k_4)$	RK4	$O(h^5)$
$h < 2/\|\lambda\|$	欧拉法稳定性条件	—

参考文献

Perko, L. (2001). Differential Equations and Dynamical Systems (3rd ed., Chapter 1). Springer.
Butcher, J. C. (2016). Numerical Methods for Ordinary Differential Equations (3rd ed.). Wiley. — Runge-Kutta方法的权威参考
Press, W. H., et al. (2007). Numerical Recipes (3rd ed., Chapter 17). Cambridge University Press.
Runge, C. (1895). “Über die numerische Auflösung von Differentialgleichungen”. — RK方法的原始论文
Kutta, W. (1901). “Beitrag zur näherungsweisen Integration totaler Differentialgleichungen”. — RK4的原始论文

下一节：偏微分方程引论：从热方程到拉普拉斯方程

从 ODE 走向 PDE：热传导方程、波动方程、拉普拉斯方程。分离变量法、特征函数展开。重点落在拉普拉斯方程——衔接空气动力学的势流理论。