亚当斯—巴什福思法

亚当斯—巴什福思法（Adams–Bashforth method）是求解常微分方程初值问题的一族显式线性多步法。它把已经算出的几个时刻的导数连成插值多项式，再对这个多项式积分，预测下一个状态。新状态不出现在待求方程右端，因此每步不必求解隐式方程；代价是需要历史数据，并受到稳定性步长的限制。

从两个已知斜率预测下一步

考虑无量纲衰减模型 $y^{\prime }=-y$、$y(0)=1$，精确解为 $e^{-t}$。记等间距时刻为 $t_n=t_0+nh$，其中 $h>0$；数值近似为 $y_n$，已知斜率为 $f_n=f(t_n,y_n)$。

显式欧拉法只用当前斜率。若还保存前一步的斜率，就能估计斜率正在怎样变化。在本例中，随着存量降低，负斜率也逐渐接近零；把刚才的下降速度一直保持不变会忽略这一趋势。

对一般方程 $y^{\prime }=f(t,y)$，精确解满足 $$y(t_{n+1})-y(t_n)={\int }_{t_n}^{t_{n+1}}f(t,y(t))\mathrm{d}t.$$ 令 $\theta =(t-t_n)/h$。通过 $\theta =0,-1$ 两个导数值的直线是 $$p(\theta )=f_n+\theta (f_n-f_{n-1}).$$ 将它从零积分到一，得到二阶、二步 Adams–Bashforth 格式，简称 AB2： $$y_{n+1}=y_n+h(\frac{3}{2}{f}_n-\frac{1}{2}{f}_{n-1}).$$ 系数中的负号来自向未来外推，并非把前一步的状态直接减去。积分的对象是导数。Adams 方法的插值积分构造和线性多步法的记号见 Driscoll 与 Braun 的教材。

示意取导数曲线 g(t)=−exp(−t)，用 t=0、0.2 两点外推到 0.4。着色部分是导数近似的积分，不是状态曲线下的面积。

更一般地，通过 $f_n,f_{n-1},\dots ,f_{n-k+1}$ 的次数不超过 $k-1$ 的插值多项式积分，给出固定步长的 $k$ 步 Adams–Bashforth 法。例如三步格式是 $$y_{n+1}=y_n+\frac{h}{12}(23f_n-16f_{n-1}+5f_{n-2}).$$ 增加历史点并不意味着任意步长都更稳定；阶数、起步误差和稳定区域要分别检查。

手算衰减过程与起步

AB2 在计算 $y_2$ 前需要 $y_0,y_1$。这里用显式梯形起步：先预测 ${\hat{y}}_1=y_0+hf(t_0,y_0)$，再平均初始斜率与预测点斜率。取 $h=0.2$，得到 $${\hat{y}}_1=0.8,y_1=1+0.1(-1-0.8)=\mathrm{0.82.}$$ 因此 $f_0=-1,f_1=-0.82$。下一步为 $$y_2=0.82+0.2[-\frac{3}{2}(0.82)+\frac{1}{2}]=\mathrm{0.674.}$$ 再用 $f_2=-0.674$ 推进： $$y_3=0.674+0.2[-\frac{3}{2}(0.674)+\frac{1}{2}(0.82)]=\mathrm{0.5538.}$$ 精确值分别为 $e^{-0.4}\approx 0.670320046$ 和 $e^{-0.6}\approx 0.548811636$。这组结果不是将误差藏在首步：起步值也来自同一个可复算算法。

一个有解的变式是 $y^{\prime }=t,y(0)=0$。显式梯形起步得到 $y_1=h^2/2$。由于真实导数本来就是直线，后续插值没有误差，AB2 给出 $y_n=t_n^2/2$。例如 $h=0.2$ 时，$y_1=0.02,y_2=0.08$。它说明二阶方法能够精确处理某些特殊问题，而非所有误差都必然非零。

二阶精度的误差来源

先把精确历史值放进格式，以排除已有误差。设沿精确解的导数函数 $g(t)=y^{\prime }(t)$ 有连续二阶导数，且 $|g^{″}|\le M$。线性插值余项在下一段满足 $$g(t_n+\theta h)-p(\theta )=\frac{h^2}{2}\theta (\theta +1)g^{″}({\zeta }_{\theta }),0\le \theta \le 1.$$ 积分并取绝对值，单步状态缺陷因而不超过 $$\frac{M{h}^3}{2}{\int }_0^1\theta (\theta +1)\mathrm{d}\theta =\frac{5}{12}M{h}^3.$$ 这里的三次方是一步状态误差；约 $1/h$ 步的累积通常留下二次方的全局误差。

这个累积可以证明。若 $f$ 对状态满足常数为 $L$ 的利普希茨条件，记状态误差为 $e_n$，上述缺陷界为 $C{h}^3$。两份递推相减得到 $$|e_{n+1}|\le (1+3Lh/2)|e_n|+(Lh/2)|e_{n-1}|+C{h}^3.$$ 令 $E_n=\underset{0\le j\le n}{max}|e_j|$，便有 $E_{n+1}\le (1+2Lh)E_n+C{h}^3$。在固定终止时间 $T$ 内迭代，并用 $(1+2Lh{)}^n\le e^{2LT}$，可得 $$E_n\le e^{2LT}(E_1+CT{h}^2).$$ 所以准确到 $O(h^2)$ 的起步值，加上足够光滑的解，保证全局二阶收敛。证明也标出了前提：历史误差会被继承，起步不能随意指定。

衰减问题中的稳定步长

对测试方程 $y^{\prime }=\lambda y$，令 $z=h\lambda$，格式成为 $$y_{n+1}=(1+3z/2)y_n-(z/2)y_{n-1}.$$ 试取模式 $y_n=r^n$，得到特征方程 $$r^2-(1+3z/2)r+z/2=0.$$ 对负实数 $z$，两个根都严格位于单位圆内的区间为 $-1<z<0$。端点 $z=-1$ 出现根 $-1$，扰动可以不断换号而不衰减。若本例用 $h=2$，同样的起步给出 $y_1=1$，后续值为 $y_2=-1,y_3=3,y_4=-7$：真实解趋零，离散值却发散。精度阶数没有排除这种不稳定。

这也是刚性微分方程需要专门积分器的原因。软件中的变步长、变阶 Adams 算法还会重新计算系数、估计误差并选择步长；不能把固定系数 AB2 直接称为完整的自适应求解器。SUNDIALS 对实际 Adams 与 BDF 求解器的区分说明了这一点。

可运行的标量实现

下列代码对标量右端使用 AB2，首步用上述显式梯形法。输入为函数、有限初值、正步长和非负整数步数；返回包含初值的时间与状态列表，不修改调用者的数据。每个常规步只需一次新的右端求值；若把一次右端求值视为常数成本，总时间为 $O(N)$，保存全部轨迹占 $O(N)$ 空间，滚动历史本身只需常数空间。

Python 3

import math


def ab2(f, t0, y0, h, steps):
    if not isinstance(steps, int) or steps < 0:
        raise ValueError("steps must be a nonnegative integer")
    t0, y0, h = map(float, (t0, y0, h))
    if not all(math.isfinite(v) for v in (t0, y0, h)) or h <= 0:
        raise ValueError("finite data and positive h required")

    def slope(t, y):
        value = float(f(t, y))
        if not math.isfinite(value):
            raise ArithmeticError("nonfinite derivative")
        return value

    ts, ys = [t0], [y0]
    if steps == 0:
        return ts, ys
    previous = slope(t0, y0)
    predicted = y0 + h * previous
    first_time = t0 + h
    if not math.isfinite(first_time) or not math.isfinite(predicted):
        raise OverflowError("nonfinite starting prediction")
    y = y0 + h * (previous + slope(first_time, predicted)) / 2
    if not math.isfinite(y):
        raise OverflowError("nonfinite state")
    ts.append(first_time)
    ys.append(y)
    if steps == 1:
        return ts, ys
    current = slope(ts[-1], y)
    for n in range(1, steps):
        y = ys[-1] + h * (1.5 * current - 0.5 * previous)
        t = t0 + (n + 1) * h
        if not math.isfinite(t) or not math.isfinite(y):
            raise OverflowError("nonfinite time or state")
        ts.append(t)
        ys.append(y)
        previous, current = current, slope(t, y)
    return ts, ys


ts, ys = ab2(lambda t, y: -y, 0, 1, 0.2, 3)
print([round(v, 6) for v in ys])  # [1.0, 0.82, 0.674, 0.5538]

有限输入仍不保证浮点运算永不溢出，也不保证所选步长稳定。这里检查非有限结果，但不实现误差容限或自动步长控制。对变步长数据，历史节点不再等距，必须重建相应插值权重。