龙格—库塔法

龙格—库塔法（Runge–Kutta methods，RK方法）在一个时间步内评价若干个阶段斜率，再把它们组合成常微分方程初值问题的下一步近似。它只需要当前解值，不像多步法那样直接使用多个过去时刻的解。阶段数和精度阶数是不同概念，增加阶段并不自动得到同样高的阶数。

设 $y^{'} = f (t, y)$ ，已知 $y_{n}$ 近似 $y (t_{n})$ 。最简单的显式欧拉法只取起点斜率 $f (t_{n}, y_{n})$ 。本篇用两阶段Heun方法展示RK的构造、阶数证明和可运行实现，再给出常见的经典四阶格式。

阶段与权重

一个 $s$ 阶段RK方法写为 $k_{i} = f (t_{n} + c_{i} h, y_{n} + h \sum_{j = 1}^{s} a_{i j} k_{j}), y_{n + 1} = y_{n} + h \sum_{i = 1}^{s} b_{i} k_{i} .$ 其中 $h$ 是步长， $a_{i j}, b_{i}, c_{i}$ 是方法的给定系数。若 $a_{i j} = 0$ 对所有 $j \geq i$ 成立，各阶段可按顺序直接计算，称显式方法；否则一般需联立求解阶段方程，属于隐式方法。本文讨论实标量初值问题，向量版本使用同一结构。

这里采用端点斜率平均的二阶Heun格式（也称显式梯形法）。不同教材的名称并不统一，下面的系数与公式确定本文所指的方法： $k_{1} = f (t_{n}, y_{n}), k_{2} = f (t_{n} + h, y_{n} + h k_{1}), y_{n + 1} = y_{n} + \frac{h}{2} (k_{1} + k_{2}) .$ 第二阶段在欧拉预测点处评价斜率，最终用两个斜率的平均值推进。它不是把当前值和预测值直接平均；平均的是每单位时间的变化率，再乘步长。^[1]

二阶精度的局部与全局证明

假设 $f$ 在真解附近有连续二阶偏导数，相关导数有界，时间区间固定。从一个精确点 $(t, y (t))$ 出发，沿解求导得到 $y^{″} = f_{t} + f_{y} f$ ，因此真解展开为 $y (t + h) = y + h f + \frac{h^{2}}{2} (f_{t} + f_{y} f) + O (h^{3}) .$ Heun第二阶段则满足 $f (t + h, y + h f) = f + h (f_{t} + f_{y} f) + O (h^{2}) .$ 代入两斜率平均公式，得到与真解相同的前面三项，所以每步从精确值出发的解值误差为 $O (h^{3})$ ，称三阶局部缺陷。

局部三阶还要结合稳定传播才能推出全局二阶。若 $f$ 对 $y$ 的Lipschitz常数为 $L$ ，两个起点相差 $e$ 时，第一斜率差至多 $L | e |$ ，预测点差至多 $(1 + h L) | e |$ ，第二斜率差至多 $L (1 + h L) | e |$ 。所以一步映射的输出差至多 $(1 + h L + h^{2} L^{2} / 2) | e | \leq (1 + K h) | e |$ 其中在 $0 < h \leq h_{0}$ 内可固定 $K = L + h_{0} L^{2} / 2$ 。令 $E_{n}$ 为解值误差，局部缺陷一致有界为 $C h^{3}$ ，即得 $E_{n + 1} \leq (1 + K h) E_{n} + C h^{3} .$ 迭代并利用 $(1 + K h)^{n} \leq e^{K n h}$ ，在 $n h \leq T$ 时有 $E_{n} \leq e^{K T} (E_{0} + C T h^{2})$ 。精确初值下 $E_{0} = 0$ ，故全局误差为 $O (h^{2})$ 。光滑性和Lipschitz控制分别承担局部展开与误差传播的作用。

指数增长的完整一步和变式

对 $y^{'} = y, y (0) = 1$ ，取 $h = 1 / 2$ 。第一阶段 $k_{1} = 1$ ；欧拉预测值为 $3 / 2$ ，第二阶段斜率也是 $3 / 2$ ；所以 $y_{1} = 1 + \frac{1 / 2}{2} (1 + 3 / 2) = \frac{13}{8} .$ 第二步同样将当前值乘 $1 + h + h^{2} / 2 = 13 / 8$ ，得到 $y_{2} = 169 / 64 = 2.640625$ ，可与真值 $e \approx 2.71828$ 比较。

图中将步长改为一，以分开显示各阶段：预测值二，两个斜率一、二，修正值 $5 / 2$ ；真值 $e$ 仍略高。图中的阶段点不是真解上的必经点。

单步从零一出发，欧拉预测一二处评价第二斜率二，Heun平均后得到一、二分之五，灰线为指数真解 — 虚线表示欧拉预测，实线表示平均斜率推进；阶段评价与最终更新承担不同作用。

变式。 若 $f (t, y) = 2$ 为常数，两次斜率均为二，Heun每步增加 $2 h$ ，从任意初值出发都精确。这个精确性来自该方程斜率不变，不表示每个光滑问题一步都能精确求解。

可运行的Heun步进

函数接受标量右端函数、初始时刻和值、固定正步长及非负步数，返回含初值的全部 $(t, y)$ 对。它不改写输入，检查函数评价是否有限；方程的光滑性并不能由这些有限次检查自动认证。

Python 3

from math import isfinite

def heun(function, t0, y0, h, steps):
    t0, y0, h = float(t0), float(y0), float(h)
    if not all(isfinite(v) for v in (t0, y0, h)) or h <= 0:
        raise ValueError("finite data and positive step required")
    if type(steps) is not int or steps < 0:
        raise ValueError("steps must be a nonnegative integer")
    y = y0
    output = [(t0, y0)]
    for n in range(steps):
        t = t0 + n*h
        k1 = float(function(t, y))
        k2 = float(function(t+h, y+h*k1))
        if not isfinite(k1) or not isfinite(k2):
            raise ValueError("nonfinite stage derivative")
        y += h*(k1+k2)/2
        if not isfinite(y):
            raise ArithmeticError("nonfinite updated state")
        output.append((t0+(n+1)*h, y))
    return output

trajectory = heun(lambda t, y: y, 0, 1, 0.5, 2)
assert trajectory == [(0.0, 1.0), (0.5, 13/8), (1.0, 169/64)]
assert heun(lambda t, y: 2, 0, 1, 0.25, 4)[-1][1] == 3
assert heun(lambda t, y: y, 0, 1, 0.1, 0) == [(0.0, 1.0)]
print(trajectory)

每步恰有两次右端评价， $N$ 步时间为 $O (N)$ ；保留全部轨迹占 $O (N)$ 空间，若仅保留当前状态则额外工作空间为 $O (1)$ 。向量问题还要计入状态维数和右端评价成本。

四阶方法与稳定性边界

经典四阶RK使用四个阶段： $\begin{matrix} k_{1} & = f (t_{n}, y_{n}), \\ k_{2} & = f (t_{n} + h / 2, y_{n} + h k_{1} / 2), \\ k_{3} & = f (t_{n} + h / 2, y_{n} + h k_{2} / 2), \\ k_{4} & = f (t_{n} + h, y_{n} + h k_{3}), \\ y_{n + 1} & = y_{n} + h (k_{1} + 2 k_{2} + 2 k_{3} + k_{4}) / 6. \end{matrix}$ 其系数满足四阶条件，在足够光滑且稳定传播的范围内全局四阶；这不是任意四阶段平均都具备的性质。对 $y^{'} = λ y$ ，可直接逐阶段代入得到稳定多项式 $1 + z + z^{2} / 2 + z^{3} / 6 + z^{4} / 24$ ，其中 $z = h λ$ ，与指数展开前五项一致。

阶数也不等于无条件稳定。Heun对同一测试方程的放大因子为 $R (z) = 1 + z + z^{2} / 2$ 。若 $λ = - 10, h = 0.3$ ，则 $z = - 3$ 、 $R = 2.5$ ，数值幅度每步增长，而真解按 $e^{- 3}$ 衰减。用于刚性方程时，显式高阶仍可能受很小稳定步长限制，通常需考虑合适的隐式方法。

参考文献

↑ Tobin A. Driscoll、Richard J. Braun：Fundamentals of Numerical Computation: Runge–Kutta methods，§6.4.3的ME2系数对应本文的端点平均格式；该书将Heun名称用于另一组二阶系数，经典四阶方法见式6.4.14及定义6.4.2。访问日期：2026-10-22。

[1] Tobin A. Driscoll、Richard J. Braun：Fundamentals of Numerical Computation: Runge–Kutta methods，§6.4.3的ME2系数对应本文的端点平均格式；该书将Heun名称用于另一组二阶系数，经典四阶方法见式6.4.14及定义6.4.2。访问日期：2026-10-22。

[1]