亚当斯—莫尔顿法

亚当斯—莫尔顿法（Adams–Moulton method）是一族包含新时刻导数的隐式线性多步法。它与亚当斯—巴什福思法都来自导数插值的积分；区别在于插值节点包含尚未算出的时刻，因此新状态通常要通过代数方程求出。隐式性可以改善稳定性，却同时引入方程求解的成本和收敛条件。

把区间两端的斜率都算进去

考虑无量纲初值问题 $y^{\prime }=-y,y(0)=1$。在时间区间 $[t_n,t_{n+1}]$ 上，精确增量等于导数的积分。若用连接两端导数的直线替代积分曲线，就得到 $$y_{n+1}=y_n+\frac{h}{2}[f(t_n,y_n)+f(t_{n+1},y_{n+1})].$$ 其中 $h=t_{n+1}-t_n>0$。这就是梯形法，也是二阶 Adams–Moulton 格式。本文把“二阶”用于精度；有些资料按历史步数编号，同一个格式会称为一阶步数的 Adams–Moulton 法，不能只比较编号而不比较公式。Adams 插值积分的构造见 Driscoll 与 Braun 的教材。

这个式子不是把未知量写在右边就停止计算。对当前衰减方程，代入右端可直接解出 $$(1+h/2)y_{n+1}=(1-h/2)y_n,y_{n+1}=\frac{1-h/2}{1+h/2}{y}_n.$$ 取 $h=0.2$，每步乘以 $9/11$，因此 $$y_1=\frac{9}{11}\approx 0.818181818,y_2=\frac{81}{121}\approx \mathrm{0.669421488.}$$ 精确值是 $e^{-0.2}\approx 0.818730753$、$e^{-0.4}\approx 0.670320046$。误差不是零，但全过程只解了一次线性关系，没有另加一个外推的起步值。

取真实导数 g(t)=−exp(−t)，在 0 到 0.4 上连接两端。与显式外推不同，右端导数的位置属于当前待求时间层。

高阶格式的权重来自插值

若再使用历史导数 $f_{n-1}$，在相对坐标 $\theta =(t-t_n)/h$ 的节点 $-1,0,1$ 上作二次插值，三个基函数分别是 $$\frac{\theta (\theta -1)}{2},1-{\theta }^2,\frac{\theta (\theta +1)}{2}.$$ 它们从零到一的积分为 $-1/12,2/3,5/12$，所以三阶格式为 $$y_{n+1}=y_n+\frac{h}{12}(5f_{n+1}+8f_n-f_{n-1}).$$ 这说明权重与节点有关，并不是从梯形法任意调整系数得到。三阶格式还需要起步历史。本文后面的代码与稳定性证明针对二阶梯形格式，不能直接移用于这一三阶格式。

精度与一步误差

沿精确解记 $g(t)=y^{\prime }(t)$，假设其二阶导数连续。端点线性插值的余项为 $$g(t_n+\theta h)-p(\theta )=\frac{h^2}{2}\theta (\theta -1)g^{″}({\zeta }_{\theta }).$$ 若 $|g^{″}|\le M$，积分后单步状态缺陷不超过 $$\frac{M{h}^3}{2}{\int }_0^1\theta (1-\theta )\mathrm{d}\theta =\frac{M{h}^3}{12}.$$ 这个结论用的是精确的两端值。实际计算还须考虑旧误差与隐式方程的求解误差。

若 $f$ 对状态有利普希茨常数 $L$，且 $hL/2<1$，两份格式相减并取绝对值，得到 $$(1-hL/2)|e_{n+1}|\le (1+hL/2)|e_n|+C{h}^3.$$ 在固定终止时间内反复使用这一式子，乘积因子至多按时间的指数增长，缺陷累计为 $O(h^2)$。这给出充分光滑、隐式方程准确求解时的全局二阶精度。条件 $hL/2<1$ 是此处误差估计与简单不动点求解的充分条件；它不是梯形法在所有问题上的必要稳定条件。

变式 $y^{\prime }=t,y(0)=0$ 可以手算验证。两端导数的平均为 $(t_n+t_{n+1})/2$，每步积分恰好等于 $(t_{n+1}^2-t_n^2)/2$，所以数值解精确为 $t_n^2/2$。这同时检查了右端时间不能仍写成旧时刻。

左半平面稳定与强衰减的差别

对测试方程 $y^{\prime }=\lambda y$，放大因子为 $$R(z)=\frac{1+z/2}{1-z/2},z=h\lambda .$$ 若 $\mathrm{Re}z<0$，则 $$|1-z/2{|}^2-|1+z/2{|}^2=-2\mathrm{Re}z>0.$$ 因此分母非零，而且 $|R(z)|<1$。这证明梯形法是 A 稳定的：线性衰减测试问题不会因为步长过大而发生指数放大。

但当 $z$ 沿负实轴趋于负无穷时，$R(z)\to -1$，并不趋于零。极快衰减的真实分量可能变成缓慢换号的离散分量。例如 $z=-100$ 时，$R=-49/51\approx -0.960784$，一步后仍保留近乎原来的大小。A 稳定并不等于没有振荡或能迅速消除刚性瞬态；这一区分在刚性微分方程中尤其重要。

解隐式方程与一次校正

对一般非线性右端，把新状态记为 $z$，求根问题是 $$G(z)=z-y_n-\frac{h}{2}[f(t_n,y_n)+f(t_{n+1},z)]=0.$$ 可用牛顿法，其分母为 $G^{\prime }(z)=1-(h/2)f_y(t_{n+1},z)$。Newton 迭代需要导数、合适初猜及非零分母；一般情形不能保证任意步长、任意初猜都收敛。

如果只先做一次显式 Euler 预测，再把这个预测值代入梯形右端校正一次，所得方法并不等于解出了隐式方程。在线性测试问题上，一次预测校正给出 $1+z+z^2/2$，恰是显式二阶格式的放大因子。取 $z=-4$，它为五；真正的隐式梯形因子却为 $-1/3$。要使用隐式格式的稳定性结论，必须确实把对应方程解到足够精度。生产求解器如何把 Adams 方法与非线性迭代结合，可参见 CVODE 的说明。

带残差检查的标量一步

代码输入标量右端 $f$、其状态导数 $f_y$、当前时间与状态、正步长；返回一个新状态，不改写输入。迭代上限防止无休止求解，残差阈值用于验收离散方程；它不自动等于真实解的误差。每步执行 $K$ 次 Newton 迭代时，右端及导数成本为 $O(K)$，工作空间为常数。

Python 3

import math


def trapezoid_step(f, df, t, y, h, tol=1e-12, max_iter=30):
    if not all(math.isfinite(v) for v in (t, y, h, tol)):
        raise ValueError("finite inputs required")
    if h <= 0 or tol <= 0 or not isinstance(max_iter, int) or max_iter < 1:
        raise ValueError("positive h, tolerance and iteration limit required")
    old = float(f(t, y))
    z = y + h * old
    for _ in range(max_iter):
        new = float(f(t + h, z))
        residual = z - y - h * (old + new) / 2
        if not all(math.isfinite(v) for v in (z, new, residual)):
            raise ArithmeticError("nonfinite iteration")
        if abs(residual) <= tol * (1 + abs(z)):
            return z
        jacobian = 1 - h * float(df(t + h, z)) / 2
        if not math.isfinite(jacobian) or abs(jacobian) < 1e-14:
            raise ArithmeticError("singular or ill-scaled Newton derivative")
        z -= residual / jacobian
    raise RuntimeError("Newton iteration did not reach the residual tolerance")


y = 1.0
for n in range(2):
    y = trapezoid_step(lambda t, y: -y, lambda t, y: -1, n * 0.2, y, 0.2)
    print(round(y, 9))  # 0.818181818, then 0.669421488

非线性方程可能有多个根；收敛到哪个根需要结合连续分支、初猜和步长判断。算法的高阶插值公式、非线性求解误差和真实初值问题的误差是三个层次，应分别控制。