Runge–Kutta 方法

形式陈述 ​

对初值问题 $y^{\prime }=f(t,y)$，一个 $s$-stage Runge–Kutta 方法在步长 $h_n$ 内计算

$$k_i=f(t_n+c_i{h}_n,y_n+h_n\sum _{j=1}^s{a}_{ij}{k}_j),i=1,\dots ,s,$$

再更新

$$y_{n+1}=y_n+h_n\sum _{i=1}^s{b}_i{k}_i.$$

系数用 Butcher tableau 记录：

$$\begin{array}{cc}c& A\\ & b^{\mathsf{T}}\end{array}.$$

通常要求内部一致性

$$c_i=\sum _{j=1}^s{a}_{ij}.$$

若 $a_{ij}=0$ 对所有 $j\ge i$ 成立，stage 可按 $i=1,\dots ,s$ 依次直接计算，方法是显式 RK。若 $A$ 为下三角且对角可非零，得到 diagonally implicit RK，每个 stage 依次解一个状态维数为 $d$ 的隐式方程；一般稠密 $A$ 则把全部 $sd$ 个 stage 耦合成一个非线性系统。stages 因而不是任意并行、彼此独立的函数采样，其依赖由 $A$ 精确规定。

阶数来自把数值一步映射与精确流的多元Taylor 展开公理库多元 Taylor 定理Multivariable Taylor theorem用各阶导数给出多元函数的局部多项式展开及余项控制。逐项匹配。在 $c=A\mathbf{1}$ 约定下，一至三阶的必要阶条件依次包括

$$\sum _i{b}_i=1,$$$$\sum _i{b}_i{c}_i=\frac{1}{2},$$

以及

$$\sum _i{b}_i{c}_i^2=\frac{1}{3},\sum _{i,j}{b}_i{a}_{ij}{c}_j=\frac{1}{6}.$$

第一式给一致性；满足到第 $p$ 阶的全部条件，才可在相应光滑性与误差传播条件下得到局部一步缺陷 $O(h^{p+1})$ 和固定终点上的全局误差 $O(h^p)$。高阶条件按不同导数组合分叉，可由 rooted trees 系统组织；基础算法不能凭 stage 数或几条低阶等式猜测更高阶。

对测试方程 $y^{\prime }=\lambda y$，令 $z=h\lambda$。消去 stages 得稳定函数

$$R(z)=1+z{b}^{\mathsf{T}}(I-zA{)}^{-1}\mathbf{1},y_{n+1}=R(z)y_n.$$

显式 RK 的 $A$ 严格下三角，$R$ 是有限次多项式；隐式 RK 通常得到有理函数。$R$ 的具体衰减区域与方法阶数不同，高阶显式方法仍只拥有有限稳定区域。

算法输入应包括 IVP、时间网格和 tableau；自适应误差控制与隐式 stage 求解再按需提供各自容差和工作预算。输出包括状态、函数与 Jacobian 求值数、非线性迭代数及退出状态。显式 $s$-stage 方法直接实现时每步至多需要 $s$ 次新的函数求值；若方法具有 FSAL 性质且前一步成功接受，可以复用端点导数并少算一次。向量组合成本通常为 $O(s^{2}d)$，稀疏 $A$ 可相应减少；隐式方法还需反复求解 $sd$ 维耦合系统或 $s$ 个顺序 stage 系统，成本取决于 Jacobian 和线性代数结构。

固定网格显式方法在所有 stages 返回有限值并完成更新后结束该步。隐式方法还必须让 stage 方程的尺度化残差达到容差；若非线性求解失败、函数值非有限、时间加法不再前进或预算耗尽，必须返回失败。嵌入式误差估计、步长接受与控制器属于自适应 RK 的后续层，不应隐含在某个 tableau 名称中。

直觉 ​

Euler 只在一步起点读取一个斜率，RK 方法则在同一步中安排多个试探位置。早期 stage 预测轨道会走向哪里，后续 stage 在这些预测位置重新读取斜率，最终按 $b_i$ 汇总。tableau 是这张依赖图的紧凑坐标，而不是一串可随意互换的系数。

阶条件要求这些试探斜率在所有低阶导数结构上共同复制精确流。两个方法可以使用相同 stage 数并具有不同阶，也可以阶数相同却拥有不同误差常数、存储需求和稳定函数；“更多 stage”不是单调的质量标签。

例子与边界 ​

显式 midpoint 的 tableau 为

$$\begin{array}{ccc}0& 0& 0\\ \frac{1}{2}& \frac{1}{2}& 0\\ & 0& 1\end{array},$$

即

$$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.$$

Heun 二阶法则取

$$\begin{array}{ccc}0& 0& 0\\ 1& 1& 0\\ & \frac{1}{2}& \frac{1}{2}\end{array}.$$

它先用 $k_1$ 预测终点斜率，再平均起终斜率。两者都有两个 stages，且满足二阶条件，但采样位置、权重和局部误差常数不同；“同阶”不表示逐步结果相同。

classical RK4 使用

$$\begin{array}{ccccc}0& 0& 0& 0& 0\\ \frac{1}{2}& \frac{1}{2}& 0& 0& 0\\ \frac{1}{2}& 0& \frac{1}{2}& 0& 0\\ 1& 0& 0& 1& 0\\ & \frac{1}{6}& \frac{1}{3}& \frac{1}{3}& \frac{1}{6}\end{array}$$

并对测试方程给出

$$R(z)=1+z+\frac{z^2}{2}+\frac{z^3}{6}+\frac{z^4}{24}.$$

这正是 $e^z$ 到四次项的截断，验证了 RK4 在线性测试方程上的四阶匹配。一般光滑非线性问题的四阶准确性仍来自前述完整的 rooted-tree 阶条件，不能由这一条稳定多项式单独推出。它仍是多项式，负实轴上的稳定区间有限；遇到很快衰减的刚性模式时，为保持稳定而被迫选择的步长可能远小于精度本身所需，增加阶数并未消除显式稳定限制。

隐式 RK 的 stages 也不能简单各自做一次固定点更新后接受。若 stage 残差没有达到与时间离散误差相容的尺度，代数求解误差会污染声明的阶数；全隐式方法还需处理 stage 之间的耦合 Jacobian。

推论与应用 ​

显式 midpoint、Heun 和 classical RK4 展示同一 stage 模型如何形成不同阶与稳定函数；Gauss、Radau 和 DIRK 则通过隐式系数追求不同的稳定或结构性质。这里列出它们的共同计算骨架，不把方法名称堆成推荐清单。

生产求解器常用共享 stages 的嵌入式公式同时给出高低阶近似，再根据差值控制步长。那一层需要误差尺度、接受/拒绝逻辑和控制器保护；本页只把 $A,b,c$、stage 方程、成本和失败接口固定下来。

参考资料

• John C. Butcher, Numerical Methods for Ordinary Differential Equations, 3rd ed., Wiley, 2016.
• Ernst Hairer, Syvert P. Nørsett, and Gerhard Wanner, Solving Ordinary Differential Equations I: Nonstiff Problems, 2nd rev. ed., Springer, 1993.
• NIST Digital Library of Mathematical Functions, §3.7: Ordinary Differential Equations.