线性多步法

线性多步法（linear multistep methods）用若干个过去网格点的解值和导数值共同推进常微分方程初值问题。保存历史能够减少每个新步所需的函数评价，但也让启动误差与过去的扰动进入递推。方法是否相容、扰动是否稳定传播以及给定步长下是否绝对稳定，是需要分别检查的性质。

本篇以二步Adams–Bashforth方法为主线，从导数外推推导公式，计算一个可手算的初值问题，再证明其二阶收敛，并用一个相容却不稳定的格式说明边界。

一般形式与启动

在等距网格 $t_n=t_0+nh$ 上，记 $f_n=f(t_n,y_n)$。一个 $k$ 步线性多步格式写为 $$\sum _{j=0}^k{\alpha }_j{y}_{n+j}=h\sum _{j=0}^k{\beta }_j{f}_{n+j},{\alpha }_k\ne 0.$$ 系数不依赖于 $n$。当 ${\beta }_k=0$ 时，新解值可直接算出，属于显式格式；当 ${\beta }_k\ne 0$ 时，新值还出现在右端函数内，一般需解一个方程，属于隐式格式。

初值问题通常只给一个起点值。二步法需要另一个起始值，常用合适的单步法产生。启动误差必须与目标全局精度相容；把第二个值随意指定，即使递推公式本身精度很高，也不能保证最终准确。这里始终讨论固定正步长，变步长需要重新计算系数。

积分导数外推得到AB2

精确解满足 $$y(t_{n+1})-y(t_n)={\int }_{t_n}^{t_{n+1}}f(t,y(t))dt.$$ 已知最近的两个斜率 $f_{n-1},f_n$，用通过它们的直线外推下一步内的导数。令 $\theta =(t-t_n)/h$，则前一个网格点对应 $\theta =-1$，当前点对应零，外推多项式为 $$p_f(\theta )=(1+\theta )f_n-\theta f_{n-1}.$$ 在 $0\le \theta \le 1$ 上积分，得到 $$y_{n+1}=y_n+h(\frac{3}{2}{f}_n-\frac{1}{2}{f}_{n-1}).$$ 这就是二步Adams–Bashforth方法，简称AB2。负权重来自把旧斜率的直线延长到新时间段，不是任意选取的修正项。每次推进只需增加一个新网格点的右端评价，已有斜率可以复用。

从局部缺陷到全局二阶

假设真解有连续四阶导数，并在固定时间区间上有界。将精确的过去解代入公式。在 $t_n=t$ 处展开 $$y^{\prime }(t-h)=y^{\prime }-h{y}^{″}+\frac{h^2}{2}{y}^{‴}-\frac{h^3}{6}{y}^{(4)}+o(h^3).$$ 因此格式右端给出的下一步值为 $$y+h{y}^{\prime }+\frac{h^2}{2}{y}^{″}-\frac{h^3}{4}{y}^{‴}+O(h^4).$$ 而真解的三次项是 $h^3{y}^{‴}/6$，所以精确解与格式预测之差为 $$d_n=\frac{5}{12}{h}^3{y}^{‴}(t_n)+O(h^4).$$ 这里用的是每步解值缺陷 $d_n=O(h^3)$；有些文献将它除以 $h$ 称为局部截断误差，采用后一约定时阶为二。

还要控制误差传播。假设 $f$ 对解变量的Lipschitz常数为 $L$，记 $e_n=y_n-y(t_n)$。数值递推减去精确解递推，取绝对值得 $$|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.$$ 逐步迭代并用 $(1+2Lh{)}^n\le e^{2Lnh}$，对 $nh\le T$ 有 $$E_n\le e^{2LT}(E_1+CT{h}^2).$$ 所以只要两个启动值误差均为 $O(h^2)$，全局误差也是 $O(h^2)$。由精确初值做一次Heun步可产生 $O(h^3)$ 的第二点误差，足够满足这里的启动要求。

指数增长的手算与变式

对 $y^{\prime }=y,y(0)=1$ 取 $h=1/2$。先用Heun启动：起点斜率一，欧拉预测值 $3/2$，故 $$y_1=1+\frac{1}{4}(1+3/2)=\frac{13}{8}.$$ 随后AB2使用 $f_0=1,f_1=13/8$，得到 $$y_2=\frac{13}{8}+\frac{1}{2}(\frac{3}{2}\frac{13}{8}-\frac{1}{2})=\frac{13}{8}+\frac{31}{32}=\frac{83}{32}=\mathrm{2.59375.}$$ 这里真解在 $t=1$ 为 $e\approx 2.71828$。第二步的导数外推直线在 $\theta =1$ 取 $9/4$，其单位区间积分为 $31/16$；乘实际步长二分之一后才是状态增量 $31/32$。

阴影是导数外推多项式的积分；状态增量还需乘步长。两个历史斜率决定本步外推。

变式。 若右端恒为二，两个斜率都等于二，AB2每步增加 $2h$；Heun启动也精确，因而整个离散轨迹精确。更一般地，沿真解的导数若是时间的一次多项式，直线外推的积分没有插值误差；还需要历史值本身准确，才能将这个精确性传递给数值轨迹。

可运行的固定步长实现

下面的标量实现用一次Heun步启动，保存全部轨迹。每个新AB2步只评价一次新点的右端，过去两点的斜率保持为滚动状态。

Python 3

from math import isfinite

def adams_bashforth_two(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")
    output = [(t0, y0)]
    if steps == 0:
        return output
    previous_f = float(function(t0, y0))
    predicted = y0+h*previous_f
    predictor_f = float(function(t0+h, predicted))
    if not all(isfinite(v) for v in (previous_f, predicted, predictor_f)):
        raise ValueError("nonfinite starting stage")
    y = y0+h*(previous_f+predictor_f)/2
    if not isfinite(y):
        raise ArithmeticError("nonfinite starting value")
    output.append((t0+h, y))
    if steps == 1:
        return output
    current_f = float(function(t0+h, y))
    if not isfinite(current_f):
        raise ValueError("nonfinite derivative")
    for n in range(1, steps):
        y += h*(1.5*current_f-0.5*previous_f)
        t = t0+(n+1)*h
        if not isfinite(y) or not isfinite(t):
            raise ArithmeticError("nonfinite updated state")
        output.append((t, y))
        if n+1 < steps:
            new_f = float(function(t, y))
            if not isfinite(new_f):
                raise ValueError("nonfinite derivative")
            previous_f, current_f = current_f, new_f
    return output

trajectory = adams_bashforth_two(lambda t, y: y, 0, 1, 0.5, 2)
assert trajectory == [(0.0, 1.0), (0.5, 13/8), (1.0, 83/32)]
assert adams_bashforth_two(lambda t, y: 2, 0, 1, 0.25, 4)[-1][1] == 3
assert adams_bashforth_two(lambda t, y: 0, 0, 7, 0.1, 5)[-1][1] == 7
assert adams_bashforth_two(lambda t, y: y, 0, 1, 0.1, 0) == [(0.0, 1.0)]
print(trajectory)

计算 $N$ 步需要 $O(N)$ 次固定成本的右端评价，存储全部输出需 $O(N)$ 空间；只滚动保存当前状态和两点斜率时，额外工作空间为 $O(1)$。本实现不支持非均匀步长或向量状态。

相容与零稳定不是同一条件

对一般多步法，定义 $$\rho (z)=\sum _{j=0}^k{\alpha }_j{z}^j,\sigma (z)=\sum _{j=0}^k{\beta }_j{z}^j.$$ 代入常数真解和一次真解，分别得到基本相容条件 $\rho (1)=0$ 与 ${\rho }^{\prime }(1)=\sigma (1)$。零稳定的根条件还要求 $\rho$ 的所有根模长不超过一，且单位圆上的根都是单根。它控制微小初始扰动在固定时间区间、步数随 $h\to 0$ 增长时的传播。AB2有 $\rho (z)=z^2-z$，根为零和一，满足该条件。

相容本身不够。考虑 $$y_{n+2}-3y_{n+1}+2y_n=-h{f}_n.$$ 其 $\rho (z)=(z-1)(z-2)$、$\sigma (z)=-1$，仍有 $\rho (1)=0$ 与 ${\rho }^{\prime }(1)=-1=\sigma (1)$，但根二违反零稳定。对最简单的 $y^{\prime }=0$，若起始误差为 $y_0=0,y_1=\epsilon$，递推的显式解是 $y_n=\epsilon (2^n-1)$；即使 $\epsilon =h^2$，在固定正时刻 $n=T/h$ 也会失控。

零稳定也不保证任意给定步长下绝对稳定。对 $y^{\prime }=\lambda y$，AB2的误差模式满足 $${\xi }^2-(1+3z/2)\xi +z/2=0,z=h\lambda .$$ 若 $\lambda =-10,h=0.2$，则 $z=-2$，根为 $-1\pm \sqrt{2}$，其中一个模长大于一。真解衰减时，数值递推仍可能出现增长的振荡分量。因此，多步法需要将阶数、启动方式、零稳定与问题所需的步长稳定性一起判断。