嵌入式龙格—库塔法

嵌入式龙格—库塔法（embedded Runge–Kutta methods）复用同一组阶段斜率，计算两个不同阶数的近似值，并用它们的差作为局部误差指标，决定是否接受当前步以及下一步多长。共享阶段使误差估计不必完整重做一次积分；差值仍是渐近指标，不是无条件的误差证明。

本篇使用两阶段的Heun—Euler对。它的高阶输出为二阶Heun值，低阶输出为一阶Euler值。对 $y^{\prime }=y,y(0)=1$，试步 $h=0.2$ 时，两阶段斜率为一和 $1.2$，Euler值为 $1.2$，Heun值为 $1.22$，差为 $0.02$。若允许指标上限为 $0.01$，这一步应拒绝并缩短步长。

共享阶段、两个输出

在当前点 $(t,y)$ 计算 $$k_1=f(t,y),k_2=f(t+h,y+h{k}_1).$$ 两套权重给出 $$y_L=y+h{k}_1,y_H=y+\frac{h}{2}(k_1+k_2),\delta =y_H-y_L=\frac{h}{2}(k_2-k_1).$$ 下标 $L,H$ 分别表示低阶与高阶。两次函数评价同时服务两个输出，额外计算差值只需少量算术。常用RK23、RK45等嵌入对采用更多阶段和更高阶权重，但这种“共享阶段、比较输出”的结构相同。

差值估计的是哪个误差

假设当前 $y$ 为精确解值，$f$ 足够光滑。令 $a=f_t+f_{y}f$，则 $$k_2=f+ha+O(h^2),y(t+h)=y+hf+\frac{h^2}{2}a+O(h^3).$$ 所以 $y_H$ 的一步缺陷为 $O(h^3)$，而 $y_L$ 缺少二次项；两者差值为 $$\delta =\frac{h^2}{2}a+O(h^3).$$ 它渐近反映的是低阶Euler输出的一步误差，而不是把高阶Heun输出的三阶缺陷精确测出来。因此用它控制并传播 $y_H$ 通常比较保守。一般高阶为 $p$、低阶为 $p-1$ 的嵌入对，其差值主阶通常为 $h^p$，控制指数应与这个阶数匹配。

若高低两输出差值的首项消失，或函数在一步内部发生未被阶段点采到的变化，差值就可能低估误差。阶段计算与真解的接近依赖光滑性，不能只因程序返回一个很小的数就跳过条件。

归一化指标和步长更新

对标量状态，使用绝对、相对容差组合 $$\mathrm{s}\mathrm{c}\mathrm{a}\mathrm{l}\mathrm{e}=\mathrm{a}\mathrm{t}\mathrm{o}\mathrm{l}+\mathrm{r}\mathrm{t}\mathrm{o}\mathrm{l}max(|y|,|y_H|),\mathrm{e}\mathrm{r}\mathrm{r}=\frac{|\delta |}{\mathrm{s}\mathrm{c}\mathrm{a}\mathrm{l}\mathrm{e}}.$$ 本文要求 $\mathrm{a}\mathrm{t}\mathrm{o}\mathrm{l}>0,\mathrm{r}\mathrm{t}\mathrm{o}\mathrm{l}\ge 0$，避免零状态导致零分母。若 $\mathrm{e}\mathrm{r}\mathrm{r}\le 1$，接受新状态；否则保留原时刻与原状态，重新试算。拒绝步只改变候选步长，不能仍把不合要求的值推进到下一时刻。

在差值约为 $C{h}^2$ 的范围内，要把归一化误差调至一，可取 $$h_{\mathrm{n}\mathrm{e}\mathrm{w}}=h\eta {\mathrm{e}\mathrm{r}\mathrm{r}}^{-1/2},$$ 其中 $0<\eta <1$ 是安全系数，通常还限制步长缩放的最大、最小倍数。本篇用 $\eta =0.9$。主例首次 $\mathrm{e}\mathrm{r}\mathrm{r}=2$，于是新步长为 $0.2\cdot 0.9/\sqrt{2}\approx 0.127279$；从同一初值重试时差值为 $h^2/2=0.0081$，低于 $0.01$。

图比较两个候选步长；自动控制器可以选择并非恰好减半的新步长。

变式。 若手动把步长改试为 $0.1$，则 $y_L=1.1,y_H=1.105$，差值 $0.005$，同样可接受。这个试步比自动选择的约 $0.1273$ 更短，显示自适应控制的目标是满足指标同时避免不必要的小步。

可运行的自适应标量实现

函数返回接受的状态列表和包含拒绝步的全部试步记录。终点不早于起点，初始步长为正，最多尝试指定次数；不改写用户数据。实现会检查非有限阶段和无法推进的浮点步长，但它使用的仍是上述启发式误差指标。

Python 3

from math import isfinite

def adaptive_heun(function, t0, y0, target, initial_step,
                  atol=1e-6, rtol=1e-3, max_trials=10000):
    t, y, target, h, atol, rtol = map(float,
        (t0, y0, target, initial_step, atol, rtol))
    if not all(isfinite(v) for v in (t, y, target, h, atol, rtol)):
        raise ValueError("finite input required")
    if target < t or h <= 0 or atol <= 0 or rtol < 0:
        raise ValueError("invalid interval, step or tolerances")
    if type(max_trials) is not int or max_trials < 1:
        raise ValueError("max_trials must be positive")
    states, trials = [(t, y)], []
    while t < target:
        if len(trials) >= max_trials:
            raise RuntimeError("trial limit reached")
        remaining = target - t
        h = min(h, remaining)
        endpoint = target if h == remaining else t+h
        if endpoint == t:
            raise ArithmeticError("step cannot advance time")
        k1 = float(function(t, y))
        k2 = float(function(endpoint, y+h*k1))
        low = y+h*k1
        high = y+h*(k1+k2)/2
        if not all(isfinite(v) for v in (k1, k2, low, high)):
            raise ArithmeticError("nonfinite trial")
        delta = abs(high-low)
        error = delta/(atol+rtol*max(abs(y), abs(high)))
        accepted = error <= 1
        trials.append({"t": t, "h": h, "delta": delta,
                       "scaled_error": error, "accepted": accepted})
        if accepted:
            t, y = endpoint, high
            states.append((t, y))
        factor = 2.0 if error == 0 else min(2.0, max(0.2, 0.9*error**-0.5))
        h *= factor
    return states, trials

states, trials = adaptive_heun(lambda t, y: y, 0, 1, 0.2, 0.2,
                               atol=0.01, rtol=0)
assert not trials[0]["accepted"]
assert abs(trials[0]["delta"] - 0.02) < 1e-14
assert trials[1]["t"] == 0
assert states[-1][0] == 0.2
assert all(row["scaled_error"] <= 1 for row in trials if row["accepted"])
assert adaptive_heun(lambda t, y: y, 0, 1, 0, 0.1)[0] == [(0.0, 1.0)]
print(states)

接受时使用的 $y_H$ 正是Heun更新，阶段正确性和二阶局部展开见龙格—库塔法；拒绝时状态未变，不变量为“所有保存的推进均通过所声明的指标”。这不变量不等于真误差已被严格界住。每个试步有两次右端评价，$N$ 次尝试耗时 $O(N)$，保存全部记录耗空间 $O(N)$。若不断拒绝，程序由试步上限或无法推进检查终止报错，而不是假定所有问题都必定成功到达终点。

零指标仍可能漏掉内部变化

取 $y^{\prime }=t(t-1),y(0)=0$，一步走到一。两个阶段都评价到零，因此 $y_L=y_H=0$，差值也为零；真解却满足 $$y(1)={\int }_0^1(t^2-t)\mathrm{d}t=-\frac{1}{6}.$$ 函数在两个采样端点之间的负变化被漏掉了。这个反例说明，误差控制还需合适的初始步长、对事件或快速变化的处理，以及与问题尺度相容的光滑性假设。相对容差很小也不会改变两阶段恰好都采到零的事实。