Runge–Kutta 方法 | Codeboy's Blog

形式陈述

作为常微分方程时间步进方法，取整数 $s\ge1$ ，一个 $s$ -stage Runge–Kutta 格式对初值问题 $y'=f(t,y)$ 在步长 $h_n$ 内计算

k_i =f\!\left( t_n+c_ih_n,\, y_n+h_n\sum_{j=1}^s a_{ij}k_j \right), \qquad i=1,\ldots,s,

再更新

y_{n+1} =y_n+h_n\sum_{i=1}^s b_ik_i.

系数用 Butcher tableau 记录：

\begin{array}{c|c} c&A\\ \hline &b^{\mathsf T} \end{array}.

通常要求内部一致性

c_i=\sum_{j=1}^s a_{ij}.

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

阶数来自把数值一步映射与精确流的多元Taylor 展开逐项匹配。在 $c=A\mathbf1$ 约定下，一至三阶的必要阶条件依次包括

\sum_i b_i=1,

\sum_i b_ic_i=\frac12,

以及

\sum_i b_ic_i^2=\frac13, \qquad \sum_{i,j}b_i a_{ij}c_j=\frac16.

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

对测试方程 $y'=\lambda y$ ，令 $z=h\lambda$ 。在 $I-zA$ 可逆时，消去 stages 得稳定函数

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

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

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

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

直觉

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

一个步长 h 内在四个试探位置取得蓝色阶段斜率 k₁ 至 k₄，并按权重合成为绿色终点 yₙ₊₁ — Runge–Kutta 的多阶段斜率

图以示意位置展示四个 stage 的试探与合成，具体取样时间由 tableau 决定；下文经典 RK4 的第二、三阶段同在步中点，第四阶段在步末。

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

例子与边界

显式 midpoint 的 tableau 为

\begin{array}{c|cc} 0&0&0\\ \tfrac12&\tfrac12&0\\ \hline &0&1 \end{array},

即

k_1=f(t_n,y_n), \qquad k_2=f\!\left(t_n+\frac h2,y_n+\frac h2k_1\right), \qquad y_{n+1}=y_n+hk_2.

Heun 二阶法则取

\begin{array}{c|cc} 0&0&0\\ 1&1&0\\ \hline &\tfrac12&\tfrac12 \end{array}.

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

classical RK4 使用

\begin{array}{c|cccc} 0&0&0&0&0\\ \tfrac12&\tfrac12&0&0&0\\ \tfrac12&0&\tfrac12&0&0\\ 1&0&0&1&0\\ \hline &\tfrac16&\tfrac13&\tfrac13&\tfrac16 \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 方程、成本和失败接口固定下来。

精度阶条件之外，RK 保辛条件要求阶段系数满足成对恒等式，隐式中点是最小例子。若要系统检查更高阶精度，B-级数与根树把不同嵌套导数拆成独立条件；已有对称方法也可通过组合积分器升阶消除领先误差。

参考资料

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, accessed 2026.

形式陈述 ​

反向链接

形式陈述