040 Multistep Methods 多步法

在前一节我们讨论了求初值问题 $$y'=f(t,y),y(t_0)=y_0\tag{1}$$ 的近似解的过程。我们使用 $t=t_n$ 时的数值计算下一个点 $t=t_{n+1}$ 时解 $\phi(t)$ 的近似值。也就是说，在计算 $\phi$ 在某点处的近似值仅仅依赖之前的一个点的值。因此这称为单步法（one-step method）。不过一旦有了 $t_0$ 开始的一些点解的近似值，我们不禁要问能不能利用这些信息来求解下一个点的近似值？也就是说如果在点 $t_1,t_2,\cdots,t_n$ 处近似值 $y_1,y_2,\cdots,y_n$ 已经得到了，如何利用这些信息来计算计算 $t_{n+1}$ 处的值 $y_{n+1}$？使用这些信息而不仅是最后一个点的信息的方法称为多步法（multistep method）。这一节会讨论两种多步法：亚当斯法和后向差分法。取决于依赖前面多少点的数据，我们能够达到任意精度。为了简单起见，这里的讨论假定步长 $h$ 是常量。

亚当斯法

之前讨论过初值问题 $(1)$ 满足 $$\phi(t_{n+1})-\phi(t_n)=\int_{t_n}^{t_{n+1}}\phi'(t)dt\tag{2}$$ 亚当斯法的核心思想是时用 $k$ 阶多项式 $P_k(t)$ 来近似 $\phi'(t)$。$P_k(t)$ 的系数由前 $k+1$ 个点的数据决定。

比如这里假定一阶多项式 $P_1(t)=At+B$，那么需要前面两个点 $(t_n,y_n),(t_{n-1},y_{n-1})$ 的数据。为了使 $P_1$ 在 $t=t_n$ 和 $t=t_{n-1}$ 处都能对 $\phi'$ 进行插值，要求 $P_1(t_n)=\phi'(t_n)=f(t_n,y_n)$ 以及 $P_1(t_{n-1})=\phi'(t_{n-1})=f(t_{n-1},y_{n-1})$ 同时成立。对于整数 $j$ 可以用 $f_j$ 表示 $f(t_j,y_j)$，因此 $A,B$ 必须满足 $$\begin{aligned} At_n+B&=f_n\\ At_{n-1}+B&=f_{n-1} \end{aligned}\tag{3}$$ 解的 $$A=\frac{f_n-f_{n-1}}{h},B=\frac{f_{n-1}t_n-f_nt_{n-1}}{h}\tag{4}$$ 用 $P_1(t)$ 替换 $\phi'(t)$ 并求解积分 $(2)$ 得到 $$\phi(t_{n+1})-\phi(t_n)=\int_{t_n}^{t_{n+1}}(At+B)dt=\frac{A}{2}(t_{n+1}^2-t_n^2)+B(t_{n+1}-t_n)$$ 使用 $y_n,y_{n+1}$ 替换 $\phi(t_n),\phi(t_{n+1})$，然后进行一些代数运算得到 $$y_{n+1}=y_n+\frac{3}{2}hf_n-\frac{1}{2}hf_{n-1}\tag{5}$$ 公式 $(5)$ 是二阶亚当斯-巴什弗思公式（second-order Adams-Bashforth formula）。通过 $y_n$ 和 $y_{n-1}$ 表示 $y_{n+1}$ 的显式公式，其局部截断误差正比于 $h^3$。

使用零阶 $P_0(t)=f_n$ 就得到一阶亚当斯-巴什弗思公式（first-order Adams-Bashforth formula），这就是欧拉公式。因为此时 $(2)$ 被积分的式子就是 $f_n$，右边就等于 $hf_n$。

使用更高阶的多项式和更多的数据，按照之前的思路就能得到更精确的亚当斯公式。比如假定 $P_3(t)$ 是三阶多项式，系数由四个点 $(t_n,y_n),(t_{n-1},y_{n-1}),(t_{n-2},y_{n-2}),(t_{n-3},y_{n-3})$ 确定。代入 $(2)$ 中替换 $\phi'(t)$，求积分，简化结果，得到四阶亚当斯-巴什弗思公式（fourth-order Adams-Bashforth formula） $$y_{n+1}=y_{n}+\frac{h}{24}(55f_n-59f_{n-1}+37f_{n-2}-9f_{n-3})\tag{6}$$ 局部截断误差正比于 $h^5$。

对亚当斯-巴什弗思公式的推导过程稍作变化，可以得到称为亚当斯-莫尔顿公式（Adams-Moulton formula）的结论。依旧看一个二阶的例子。令一阶多项式 $Q_1(t)=\alpha t+\beta$，但是由点 $(t_n,y_n),(t_{n+1},y_{n+1})$ 确定，那么 $\alpha,\beta$ 满足 $$\begin{aligned} \alpha t_n+\beta&=f_n\\ \alpha t_{n+1}+\beta&=f_{n+1} \end{aligned}\tag{7}$$ 那么 $$\alpha=\frac{f_{n+1}-f_n}{h},\beta=\frac{f_nt_{n+1}-f_{n+1}t_n}{h}\tag{8}$$ 代入 $(2)$ 简化得到 $$y_{n+1}=y_n+\frac{1}{2}hf_n+\frac{1}{2}f(t_{n+1},y_{n+1})\tag{9}$$ 这就是二阶亚当斯-莫尔顿公式（second-order Adams-Moulton formula）。这里写作 $f(t_{n+1},y_{n+1})$ 来强调亚当斯-莫尔顿公式是隐式的，因为未知项 $y_{n+1}$ 出现在了公式两边。局部截断误差正比于 $h^3$。

一阶亚当斯-莫尔顿公式（first-order Adams-Moulton formula）就是后向欧拉公式。和之前类似，$(2)$ 的被积分式子就是 $f_{n+1}$，很容易推导出这个结论。

和之前一样，使用更高阶的多项式可以得到更精确的公式。四阶亚当斯-莫尔顿公式（fourth-order Adams-Moulton formula）的局部截断误差正比于 $h^5$ $$y_{n+1}=y_n+\frac{h}{24}(9f_{n+1}+19f_n-5f_{n-1}+f_{n-2})\tag{10}$$ 显然这也是一个隐式公式。

尽管同阶亚当斯-巴什弗思公式和亚当斯-莫尔顿公式的局部截断误差正比于同阶 $h$ 的幂次，但是亚当斯-莫尔顿公式更精确。比如 $(6),(10)$ 的局部截断误差都正比于 $h^5$，但是亚当斯-莫尔顿公式的常量系数亚当斯-巴什弗思公式的常量系数的十分之一还要小。

那么是使用更快的显式亚当斯-巴什弗思公式还是慢一点但是更精确的隐式亚当斯-莫尔顿公式呢？这取决于是否可以通过使用更精确的公式以增加步长，进而减少步数来抵消每一步的附加开销。

事实上数值分析专家试图将这两类公式结合起来，推导处了预测-修正法（predictor-corrector method），从而兼顾计算的简便性与精度。一旦知道了 $y_{n-3},y_{n-2},y_{n-1},y_n$，那么可以计算 $f_{n-3},f_{n-2},f_{n-1},y_n$，使用亚当斯-巴什弗思公式 $(6)$ 计算得到 $y_{n+1}$，那么可以得到 $f_{n+1}$，此时亚当斯-莫尔顿公式就是显式公式，得到更好的 $y_{n+1}$。如果 $y_{n+1}$ 的变化很大，那么可以继续使用 $(10)$ 进一步修正 $y_{n+1}$。不过如果需要超过一次甚至两次修正，说明步长 $h$ 太大需要减少一些。

为了使用多步法，需要计算得到一些 $y_j$。比如四阶亚当斯-莫尔顿公式需要 $y_1,y_2$，四阶亚当斯-巴什弗思公式还需要 $y_3$。一种方法是使用精度相当的单步法计算这些值。比如使用四阶多步法，可以使用四阶龙格-库塔法计算开始的几个值。下面的例 1 就是这样做的。

另一种方式是使用低阶方法配合更小的 $h$ 来计算 $y_1$，逐步增加步长和阶数直到计算得到了足够的初始值。

例 1 回到之前的初值问题 $$y'=1-t+4y,\quad y(0)=1\tag{11}$$ 使用步长 $h=0.1$，分别使用四阶亚当斯-巴什弗思公式、四阶亚当斯-莫尔顿公式和四阶预测-修正法计算 $t=0.4$ 时解 $y=\phi(t)$ 的近似值。

解：首先使用四阶龙格-库塔公式计算 $y_1,y_2,y_3$，这个在上一节的例子中已经计算得到了。这样就得到了相应的 $f(t,y)$ $$\begin{aligned} &y_0=1&&f_0=5\\ &y_1=1.6089333&&f_1=7.3357332\\ &y_2=2.5050062&&f_2=10.820025\\ &y_3=3.8294145&&f_3=16.017658\\ \end{aligned}$$ 使用四阶亚当斯-巴什弗思公式 $(6)$ 可以计算得到 $y_4=5.7836305$。相比 $t=0.4$ 时的精确值 5.7942260，误差是 -0.0105955。

四阶亚当斯-莫尔顿公式 $(10)$ 得到 $$y_4=4.9251275+0.15y_4$$ 那么 $y_4=5.7942676$，误差仅有 0.0000416。

最后，使用四阶亚当斯-巴什弗思公式的结果作为预测 $y_4=5.7836305$，那么 $f_4=23.734522$，代入 $(10)$ 修正得到 $y_4=5.7926721$，误差是 -0.0015539。

亚当斯-巴什弗思法最简单也最快，仅涉及显式公式，也最不精确。使用亚当斯-莫尔顿公式作为修正，公式仍旧是显式的，不过增加了一些计算量。这种方法将误差修正为预测值的七分之一。亚当斯-莫尔顿法最精确，比修正值误差的四十分之一。

不过亚当斯-莫尔顿法是隐式的，这意味着每一步需要求解方程。这个问题是线性问题，方程容易求解，如果是非线性的，那么这个过程会更耗时。

使用步长为 $h=0.1$ 的龙格-库塔法求解的 $y_4=5.7927853$，误差是 -0.0014407。对于这个问题，它的精度和预测-修正法的精度差不多。

后向差分公式

另一个多步法是使用多项式 $P_k(t)$ 近似初值问题 $(1)$ 中的 $\phi(t)$ 而不是它的导数 $\phi'(t)$。然后对 $P_k(t)$ 求导，令 $P_k'(t_{n+1})=f(t_{n+1},y_{n+1})$ 得到 $y_{n+1}$ 的隐式公式。这称为后向差分公式（backward differentiation formula）。

最简单的例子是一阶多项式 $P_1(t)=At+B$，那么要满足 $P_1(t_n)=y_n,P_1(t_{n+1})=y_{n+1}$，那么 $A,B$ 要满足 $$\begin{aligned} At_n+B&=y_n\\ At_{n+1}+B&=y_{n+1} \end{aligned}\tag{12}$$ 求解得到 $$A=\frac{y_{n+1}-y_n}{h},B=\frac{y_{n+1}t_n-y_nt_{n+1}}{h}\tag{13}$$ 由于 $P_1'(t)=A$，又因为 $P_1'(t_{n+1})=f(t_{n+1},y_{n+1})$，因此 $A=f(t_{n+1},y_{n+1})$。结合 $(13)$ 就得到了一阶后向差分公式（first-order backward differentiation formula） $$y_{n+1}=y_n+hf(t_{n+1},y_{n+1})\tag{14}$$ 这就是在 8.1 小节讨论的后向欧拉公式。

使用更高阶的多项式和相应更多的点，可以得到任意阶后向差分公式。下面是二阶后向差分公式（second-order backward differentiation formula） $$y_{n+1}=\frac{1}{3}(4y_n-y_{n-1}+2hf(t_{n+1},y_{n+1}))\tag{15}$$ 下面是四阶后向差分公式（fourth-order backward differentiation formula） $$y_{n+1}=\frac{1}{25}(48y_n-36y_{n-1}+16y_{n-2}-3y_{n-3}+12hf(t_{n+1},y_{n+1}))\tag{16}$$ 这两个公式的局部截断误差分别正比于 $h^3,h^5$。

例 2 使用步长 $h=0.1$ 的四阶后向差分公式来估算例 1 中在 $t=0.4$ 时解 $y=\phi(t)$ 的近似值。

解：$y_0,\cdots,y_3$ 在例 1 中已经给出，$n=3,h=0.1$，都代入 $(16)$ 得到 $$y_4=4.6837842+0.192y_4$$ 因此 $$y_4=5.7967626$$ 和精确值 $\phi(t)=5.7942260$ 相比，误差是 0.0025366。这比四阶亚当斯-巴什弗思的结果要好，但是不如四阶预测-修正法，远不如亚当斯-莫尔顿法。

比较单步法与多步法时，要综合考虑多个因素。四阶龙格-库塔法在每一步都需要求四次 $f$ 函数的值；而四阶亚当斯-巴什弗思法（除了开始阶段）每步仅需一次求值，预测-修正法需要两次。因此，在给定步长 $h$ 的情况下，后两种方法的运行速度很可能比龙格-库塔法快得多。不过，如果龙格-库塔法的精度更高，并因此可以采用更大的步长，更少的步数，那么这种速度上的差异就会缩小，甚至可能会消失。

亚当斯-莫尔顿和后向差分公式还考虑到在每一步中求解隐式方程的难度。所有多步法都存在一个潜在的缺点，开始步骤中的误差可能会反馈到后续计算中，从而产生不利结果。另一方面，多步法底层的多项式近似特性使得在需要时可以很容易地估算点之间任何位置的解。多步法变得流行的原因是它能相对容易地估算每一步的误差并调整阶数或步长来对其进行控制。