后向微分公式

后向微分公式（backward differentiation formula，BDF）是一族隐式线性多步法。它先用当前待求状态与若干历史状态作插值，再令插值多项式在新时刻的导数等于微分方程右端。它与 Adams 方法的区别是：BDF 近似的是状态的导数，Adams 方法近似的是导数的积分。

低阶 BDF 对快速衰减分量具有良好稳定性，是刚性微分方程求解器中的常用选择。不过“隐式”并不自动提供任意高阶、任意步长的稳定保证，非线性方程也仍要实际求解。

从三点插值推导 BDF2

考虑 $y^{\prime }=f(t,y)$，时间网格等距，步长为 $h>0$。先看新时刻 $t_{n+1}$ 及其前两个点，把相对坐标写成 $\theta =(t-t_{n+1})/h$。三个节点为 $-2,-1,0$。

通过状态值 $y_{n-1},y_n,y_{n+1}$ 的二次插值多项式为 $$p(\theta )=\frac{\theta (\theta +1)}{2}{y}_{n-1}-\theta (\theta +2)y_n+\frac{(\theta +1)(\theta +2)}{2}{y}_{n+1}.$$ 在 $\theta =0$ 求导，并注意 $\mathrm{d}\theta /\mathrm{d}t=1/h$： $$p_t(t_{n+1})=\frac{y_{n-1}-4y_n+3y_{n+1}}{2h}.$$ 令这个导数等于方程在新时刻的导数，得到二阶 BDF： $$\frac{3y_{n+1}-4y_n+y_{n-1}}{2h}=f(t_{n+1},y_{n+1}).$$ 一阶 BDF 则是后向 Euler 法 $(y_{n+1}-y_n)/h=f(t_{n+1},y_{n+1})$。BDF 的插值微分构造与多步法稳定性见 Driscoll 与 Braun 的教材。

示意取 y(t)=exp(−t)，三个节点为 0、0.2、0.4。BDF 使用右端切线；历史状态本身参与插值，而不是先对三个导数作平均。

衰减算例与首步误差

取 $y^{\prime }=-y,y(0)=1,h=0.2$。用后向 Euler 起步，有 $$y_1=\frac{1}{1+h}=\frac{5}{6}.$$ 从第二步起，BDF2 可解为 $$y_{n+1}=\frac{4y_n-y_{n-1}}{3+2h}.$$ 因此 $$y_2=\frac{4(5/6)-1}{3.4}=\frac{35}{51}\approx \mathrm{0.686274510.}$$ 这个结果高于精确值 $e^{-0.4}\approx 0.670320046$，其中既有 BDF2 误差，也有首步偏差。如果改用梯形首步 $y_1=9/11$，则 $$y_2=\frac{4(9/11)-1}{3.4}=\frac{125}{187}\approx \mathrm{0.668449198.}$$ 较准确的起步改善了这一次计算，但两种起步不能在同一张误差表中混为一个算法。

后向 Euler 若重复执行，整体只有一阶；但只用它起步一次，首步状态误差为 $O(h^2)$。这与二阶多步法所需的起步精度相容，因此不意味着整个 BDF2 程序必然退化为一阶。高精度或变阶软件会更细致地处理起步和误差估计。SUNDIALS 的 CVODE 使用变阶 BDF，并明确给出其实际支持的阶数范围。

一致性与历史误差

把精确解在新时刻展开： $$y(t_{n+1}-h)=y-h{y}^{\prime }+\frac{h^2}{2}{y}^{″}-\frac{h^3}{6}{y}^{‴}+O(h^4),$$ $$y(t_{n+1}-2h)=y-2h{y}^{\prime }+2h^2{y}^{″}-\frac{4h^3}{3}{y}^{‴}+O(h^4).$$ 代入差分导数，常数项与二阶导数项抵消，留下 $$\frac{3y-4y(t_{n+1}-h)+y(t_{n+1}-2h)}{2h}=y^{\prime }-\frac{h^2}{3}{y}^{‴}+O(h^3).$$ 所以导数残差为 $O(h^2)$，乘以一个步长后对应的状态缺陷为 $O(h^3)$。这里需要相应时间导数存在并有界；不光滑输入会破坏直接使用这个展开的条件。

还要检查历史误差会不会自行放大。令右端为零，递推的特征多项式是 $$3r^2-4r+1=(r-1)(3r-1).$$ 于是任意历史扰动的齐次解为 $e_n=C_1+C_2(1/3{)}^n$。常量误差会保留，但第二个寄生模式衰减，不会仅因步数增加而无界增长。这是零稳定性的具体表现。结合一致性、有限区间上的利普希茨条件及合适起步，可得二阶收敛。

一个直接可检验的变式是 $y^{\prime }=t$。若使用精确起步 $y_0=0,y_1=h^2/2$，由于二次状态 $t^2/2$ 的插值就是它自身，BDF2 从此精确再现这个解。若首步故意偏移，前述常量与衰减模式会解释后续的偏差，不能把偏差全归于导数公式。

稳定性的平方恒等式

对衰减方程 $y^{\prime }=-\mu y$、$\mu \ge 0$，BDF2 满足 $3a-4b+c=-2h\mu a$，其中 $a=y_{n+1},b=y_n,c=y_{n-1}$。直接展开平方可验证 $$2a(3a-4b+c)=a^2+(2a-b{)}^2-b^2-(2b-c{)}^2+(a-2b+c{)}^2.$$ 定义包含两层历史的量 $E_n=y_n^2+(2y_n-y_{n-1}{)}^2$，就得到 $$E_{n+1}-E_n=-4h\mu y_{n+1}^2-(y_{n+1}-2y_n+y_{n-1}{)}^2\le 0.$$ 它证明任意正步长都不会把这种线性衰减递推变成无界增长。这里控制的是两层状态组成的二次量，并非断言每个单独状态都保持正值或逐步减小。

对于复测试参数 $\lambda$，把乘积改成复内积的实部，平方改成模平方，同一恒等式给出右端 $4h\mathrm{Re}\lambda |y_{n+1}{|}^2$。因此 BDF2 在左半平面 A 稳定。这个结论属于 BDF1、BDF2；不能从低阶例子推出所有高阶 BDF 都 A 稳定。常步长 BDF 的三至六阶具有不同的稳定扇区，更高阶还会失去零稳定性。

线性衰减族的可运行实现

代码只求 $y^{\prime }=-\mu y$，以一个后向 Euler 步起步。输入为非负 $\mu$、有限初值、正步长和非负整数步数；返回包含初值的列表，不修改输入。每步是常数次运算，总时间 $O(N)$，保存结果占 $O(N)$，两层历史只需常数空间。对一般非线性 $f$，须另外实现并检查新时刻方程的求根。

Python 3

import math


def bdf2_decay(mu, y0, h, steps):
    if not isinstance(steps, int) or steps < 0:
        raise ValueError("steps must be a nonnegative integer")
    mu, y0, h = map(float, (mu, y0, h))
    if not all(math.isfinite(v) for v in (mu, y0, h)) or mu < 0 or h <= 0:
        raise ValueError("finite data, mu nonnegative and h positive required")
    z = mu * h
    if not math.isfinite(z) or not math.isfinite(3 + 2 * z):
        raise OverflowError("step coefficient outside floating-point range")
    values = [y0]
    if steps == 0:
        return values
    values.append(y0 / (1 + z))
    for _ in range(1, steps):
        new = (4 * values[-1] - values[-2]) / (3 + 2 * z)
        if not math.isfinite(new):
            raise OverflowError("intermediate or result outside floating-point range")
        values.append(new)
    return values


print(bdf2_decay(1, 1, 0.2, 2))  # [1.0, 0.8333333333333334, 0.6862745098039217]
print(bdf2_decay(0, 3, 0.2, 3))  # [3.0, 3.0, 3.0, 3.0]

常步长系数来自等距节点。自适应步长改变了节点间距，不能原样沿用 $3,-4,1$；频繁或过大的步长比变化也会影响多步法的稳定性。实际计算还需要误差估计、容限、起步和非线性求解控制。