线性多步法

线性多步法是一种初值问题的常微分方程数值方法。在概念上，初值问题的数位方法是由初值开始，在时间上前进一小段，求解下一个点的数值。重复此步骤，直到要求解的区间都完成为止。

单步法（像欧拉方法）只用前一个点以及其微分来计算目前的值，单步法的龙格-库塔法会用前一个点以及和中间点中间的值来计算目前的值，不过在计算下一个值时，就不考虑之前的资讯。

多步法和单步法的差异是，多步法会考虑以前的资料，以提升效率。因此。多步法会参考之前的数个点以及数点微分值。线性多步法会使用之前的点和微分值的线性组合。

定义

常微分方程的数值方法设法找到以下初值问题的近似解 $${\displaystyle y'=f(t,y),\quad y(t_{0})=y_{0}.}$$

结果是在各离散时间${\displaystyle t_{i}}$下，${\displaystyle y(t)}$的近似值： $${\displaystyle y_{i}\approx y(t_{i})\quad {\text{where}}\quad t_{i}=t_{0}+ih,}$$ 其中${\displaystyle h}$是时间步长（有时会表示为${\displaystyle \Delta t}$），而${\displaystyle i}$是整数。

多步法会用以前${\displaystyle s}$步以内的资讯来计算下一个值。而线性多步法会用${\displaystyle y_{i}}$和${\displaystyle f(t_{i},y_{i})}$的线性组合，来计算${\displaystyle y}$。因此，线性多步法可以表示如下 $${\displaystyle {\begin{aligned}&y_{n+s}+a_{s-1}\cdot y_{n+s-1}+a_{s-2}\cdot y_{n+s-2}+\cdots +a_{0}\cdot y_{n}\\&\qquad {}=h\cdot \left(b_{s}\cdot f(t_{n+s},y_{n+s})+b_{s-1}\cdot f(t_{n+s-1},y_{n+s-1})+\cdots +b_{0}\cdot f(t_{n},y_{n})\right)\\&\Leftrightarrow \sum _{j=0}^{s}a_{j}y_{n+j}=h\sum _{j=0}^{s}b_{j}f(t_{n+j},y_{n+j}),\end{aligned}}}$$ 而${\displaystyle a_{s}=1}$。系数${\displaystyle a_{0},\dotsc ,a_{s-1}}$和${\displaystyle b_{0},\dotsc ,b_{s}}$决定所使用的方法。设计者一方面要找到精确解的理想近似解，另一方面也要设计方便计算的方法。为了简化计算，中许多系数会是0。

可以用系数来区分显式和隐式方法。若${\displaystyle b_{s}=0}$，此方法则为显式法，因为可以直接计算${\displaystyle y_{n+s}}$。若${\displaystyle b_{s}\neq 0}$，则${\displaystyle y_{n+s}}$的值和${\displaystyle f(t_{n+s},y_{n+s})}$有关，需设法求解方程才能得到${\displaystyle y_{n+s}}$, 此方法为隐式法。求解方程时，常会使用迭代法（例如牛顿法）求解。

有时，显式的多步法会用来“预估”${\displaystyle y_{n+s}}$的值，之后再用隐式公式中来得到“校正”后的值，这称为预估-校正方法。

例子

考虑以下问题 $${\displaystyle y'=f(t,y)=y,\quad y(0)=1.}$$ 其精确解是${\displaystyle y(t)=e^{t}}$。

单步欧拉法

欧拉法是一种简易的单步法： $${\displaystyle y_{n+1}=y_{n}+hf(t_{n},y_{n}).}$$ 可以将欧拉法视为是退化的显式多步法，其步数只有一步。

用此方法，配合步长${\displaystyle h={\tfrac {1}{2}}}$求解问题${\displaystyle y'=y}$，得到以下的解： $${\displaystyle {\begin{aligned}y_{1}&=y_{0}+hf(t_{0},y_{0})=1+{\tfrac {1}{2}}\cdot 1=1.5,\\y_{2}&=y_{1}+hf(t_{1},y_{1})=1.5+{\tfrac {1}{2}}\cdot 1.5=2.25,\\y_{3}&=y_{2}+hf(t_{2},y_{2})=2.25+{\tfrac {1}{2}}\cdot 2.25=3.375,\\y_{4}&=y_{3}+hf(t_{3},y_{3})=3.375+{\tfrac {1}{2}}\cdot 3.375=5.0625.\end{aligned}}}$$

二步Adams–Bashforth法

二步的Adams–Bashforth法是简单的多步法 $${\displaystyle y_{n+2}=y_{n+1}+{\tfrac {3}{2}}hf(t_{n+1},y_{n+1})-{\tfrac {1}{2}}hf(t_{n},y_{n}).}$$ 此方法会用到二个值，${\displaystyle y_{n+1}}$和${\displaystyle y_{n}}$来计算下一个值${\displaystyle y_{n+2}}$。不过，初始问题只会有一个值${\displaystyle y_{0}=1}$。一种可能作法是用欧拉法计算${\displaystyle y_{1}}$，当作第二个值。在此选择下，Adams–Bashforth法得到结果如下（四舍五入到小数第四位）： $${\displaystyle {\begin{aligned}y_{2}&=y_{1}+{\tfrac {3}{2}}hf(t_{1},y_{1})-{\tfrac {1}{2}}hf(t_{0},y_{0})=1.5+{\tfrac {3}{2}}\cdot {\tfrac {1}{2}}\cdot 1.5-{\tfrac {1}{2}}\cdot {\tfrac {1}{2}}\cdot 1=2.375,\\y_{3}&=y_{2}+{\tfrac {3}{2}}hf(t_{2},y_{2})-{\tfrac {1}{2}}hf(t_{1},y_{1})=2.375+{\tfrac {3}{2}}\cdot {\tfrac {1}{2}}\cdot 2.375-{\tfrac {1}{2}}\cdot {\tfrac {1}{2}}\cdot 1.5=3.7812,\\y_{4}&=y_{3}+{\tfrac {3}{2}}hf(t_{3},y_{3})-{\tfrac {1}{2}}hf(t_{2},y_{2})=3.7812+{\tfrac {3}{2}}\cdot {\tfrac {1}{2}}\cdot 3.7812-{\tfrac {1}{2}}\cdot {\tfrac {1}{2}}\cdot 2.375=6.0234.\end{aligned}}}$$ 在${\displaystyle t=t_{4}=2}$的解是${\displaystyle e^{2}=7.3891\ldots }$，二步的Adams–Bashforth法比欧拉法准确。一般而言，只要步长够小，二步的Adams–Bashforth法都会比欧拉法准确。

多步法分类

常见的多步法可以分为三类：Adams–Bashforth法、Adams–Moulton法及反向微分公式（英语：backward differentiation formula）（BDF）。

Adams–Bashforth法

Adams–Bashforth法是显式法，系数是${\displaystyle a_{s-1}=-1}$，${\displaystyle a_{s-2}=\cdots =a_{0}=0}$，而选择${\displaystyle b_{j}}$使此方法可以到s阶（因此在同一阶下，此方法的系数唯一）。

s = 1, 2, 3, 4, 5的Adams–Bashforth法如下，其中的第一式就是前向欧拉法：(Hairer, Nørsett & Wanner 1993，§III.1; Butcher 2003，第103页)： $${\displaystyle {\begin{aligned}y_{n+1}&=y_{n}+hf(t_{n},y_{n}),\\y_{n+2}&=y_{n+1}+h\left({\frac {3}{2}}f(t_{n+1},y_{n+1})-{\frac {1}{2}}f(t_{n},y_{n})\right),\\y_{n+3}&=y_{n+2}+h\left({\frac {23}{12}}f(t_{n+2},y_{n+2})-{\frac {16}{12}}f(t_{n+1},y_{n+1})+{\frac {5}{12}}f(t_{n},y_{n})\right),\\y_{n+4}&=y_{n+3}+h\left({\frac {55}{24}}f(t_{n+3},y_{n+3})-{\frac {59}{24}}f(t_{n+2},y_{n+2})+{\frac {37}{24}}f(t_{n+1},y_{n+1})-{\frac {9}{24}}f(t_{n},y_{n})\right),\\y_{n+5}&=y_{n+4}+h\left({\frac {1901}{720}}f(t_{n+4},y_{n+4})-{\frac {2774}{720}}f(t_{n+3},y_{n+3})+{\frac {2616}{720}}f(t_{n+2},y_{n+2})-{\frac {1274}{720}}f(t_{n+1},y_{n+1})+{\frac {251}{720}}f(t_{n},y_{n})\right).\end{aligned}}}$$

系数可以用以下方法计算。用多项式插值找到${\displaystyle s-1}$阶的多项式p使得 $${\displaystyle p(t_{n+i})=f(t_{n+i},y_{n+i}),\qquad {\text{for }}i=0,\ldots ,s-1.}$$ 多项式内插的拉格朗日插值法可以得到 $${\displaystyle p(t)=\sum _{j=0}^{s-1}{\frac {(-1)^{s-j-1}f(t_{n+j},y_{n+j})}{j!(s-j-1)!h^{s-1}}}\prod _{i=0 \atop i\neq j}^{s-1}(t-t_{n+i}).}$$ 多项式p是要求解的微分方程${\displaystyle y'=f(t,y)}$等号右边的局部好近似值，因此改为考虑方程${\displaystyle y'=p(t)}$，此方法有精确解，就是p的积分。因此得到下式 $${\displaystyle y_{n+s}=y_{n+s-1}+\int _{t_{n+s-1}}^{t_{n+s}}p(t)\,\mathrm {d} t.}$$ 当上述公式的p取代之后，就是Adams–Bashforth法。其系数${\displaystyle b_{j}}$可以由下式求得 $${\displaystyle b_{s-j-1}={\frac {(-1)^{j}}{j!(s-j-1)!}}\int _{0}^{1}\prod _{i=0 \atop i\neq j}^{s-1}(u+i)\,\mathrm {d} u,\qquad {\text{for }}j=0,\ldots ,s-1.}$$ 将${\displaystyle f(t,y)}$用内插的p取代，引进了数量级为h^s的误差，而s步的Adams–Bashforth法确实是s阶(Iserles 1996，§2.1)。

Adams–Bashforth法是由约翰·柯西·亚当斯求解毛细现象建模的微分方程而来，此微分方程出自Francis Bashforth（英语：Francis Bashforth）。Bashforth (1883)发表了他的理论以及亚当斯的数值方法(Goldstine 1977)。

Adams–Moulton法

Adams–Moulton法和Adams–Bashforth法类似，也是${\displaystyle a_{s-1}=-1}$，${\displaystyle a_{s-2}=\cdots =a_{0}=0}$，也是选择系数b到最高可能阶数。但Adams–Moulton法是隐式法。去除了${\displaystyle b_{s}=0}$的限制，s步的Adams–Moulton法可以到${\displaystyle s+1}$阶，而s步的Adams–Bashforth法只有s阶。

s = 0, 1, 2, 3, 4的Adams–Moulton法如下(Hairer, Nørsett & Wanner 1993，§III.1; Quarteroni, Sacco & Saleri 2000)，其中前两个分别是反向欧拉法和梯形法则： respectively: $${\displaystyle {\begin{aligned}y_{n}&=y_{n-1}+hf(t_{n},y_{n}),\\y_{n+1}&=y_{n}+{\frac {1}{2}}h\left(f(t_{n+1},y_{n+1})+f(t_{n},y_{n})\right),\\y_{n+2}&=y_{n+1}+h\left({\frac {5}{12}}f(t_{n+2},y_{n+2})+{\frac {8}{12}}f(t_{n+1},y_{n+1})-{\frac {1}{12}}f(t_{n},y_{n})\right),\\y_{n+3}&=y_{n+2}+h\left({\frac {9}{24}}f(t_{n+3},y_{n+3})+{\frac {19}{24}}f(t_{n+2},y_{n+2})-{\frac {5}{24}}f(t_{n+1},y_{n+1})+{\frac {1}{24}}f(t_{n},y_{n})\right),\\y_{n+4}&=y_{n+3}+h\left({\frac {251}{720}}f(t_{n+4},y_{n+4})+{\frac {646}{720}}f(t_{n+3},y_{n+3})-{\frac {264}{720}}f(t_{n+2},y_{n+2})+{\frac {106}{720}}f(t_{n+1},y_{n+1})-{\frac {19}{720}}f(t_{n},y_{n})\right).\end{aligned}}}$$

Adams–Moulton法的推导类似Adams–Bashforth法，不过，用的内插法不只用${\displaystyle t_{n-1},\dots ,t_{n-s}}$的点，也会用${\displaystyle t_{n}}$点，其系数为 $${\displaystyle b_{s-j}={\frac {(-1)^{j}}{j!(s-j)!}}\int _{0}^{1}\prod _{i=0 \atop i\neq j}^{s}(u+i-1)\,\mathrm {d} u,\qquad {\text{for }}j=0,\ldots ,s.}$$

Adams–Moulton法是由约翰·柯西·亚当斯独立完成，类似Adams–Bashforth法。福里斯特·雷·莫尔顿的名字出现在方法，是因为他发现此方法可以和Adams–Bashforth法一起使用，成为一对预估-校正方法(Moulton 1926)。Milne (1926)也有类似的想法，亚当斯是用牛顿法求解隐式方程(Hairer, Nørsett & Wanner 1993，§III.1)。

反向微分公式（BDF）

反向微分公式是${\displaystyle b_{s-1}=\cdots =b_{0}=0}$的隐式法，选择其他系数使阶数可以到s阶（最大可能阶数）。此方法特别用来求得刚性方程。

分析

微分方程数值方法分析的中心概念，包括有收敛、阶数和稳定性。

收敛、阶数

第一个问题是此方法是否具有一致性：以下差分方程 $${\displaystyle {\begin{aligned}&a_{s}y_{n+s}+a_{s-1}y_{n+s-1}+a_{s-2}y_{n+s-2}+\cdots +a_{0}y_{n}\\&\qquad {}=h{\bigl (}b_{s}f(t_{n+s},y_{n+s})+b_{s-1}f(t_{n+s-1},y_{n+s-1})+\cdots +b_{0}f(t_{n},y_{n}){\bigr )},\end{aligned}}}$$ 是否是${\displaystyle y'=f(t,y)}$微分方程的良好近似？更准确的说法，多步法具有一致性，当其步阶h趋近于零时，局部截尾误差趋近于零的速度比h还快，其中局部截尾误差定义为此方法的结果${\displaystyle y_{n+s}}$，和时间${\displaystyle t_{n+s}}$时方程精确解的误差（假设所有之前的值${\displaystyle y_{n+s-1},\ldots ,y_{n}}$都是正确的）。利用泰勒级数的计算可得线性多步法具有一致性，若且唯若 $${\displaystyle \sum _{k=0}^{s-1}a_{k}=-1\quad {\text{and}}\quad \sum _{k=0}^{s}b_{k}=s+\sum _{k=0}^{s-1}ka_{k}.}$$

上述所有的方式都具有一致性(Hairer, Nørsett & Wanner 1993，§III.2)。

若方法具有一致性，下一个问题是：近似微分方程的差分方程，可以近似到多好？多步法其阶数为p若在其步阶h趋近于零时，其局部误差约为${\displaystyle O(h^{p+1})}$。这和以下方法系数的条件等效： $${\displaystyle \sum _{k=0}^{s-1}a_{k}=-1\quad {\text{and}}\quad q\sum _{k=0}^{s}k^{q-1}b_{k}=s^{q}+\sum _{k=0}^{s-1}k^{q}a_{k}{\text{ for }}q=1,\ldots ,p.}$$ s步Adams–Bashforth方法的阶数为s，而s步Adams–Moulton方法的阶数为${\displaystyle s+1}$(Hairer, Nørsett & Wanner 1993，§III.2)。

这些条件可以用特征多项式表示： $${\displaystyle \rho (z)=z^{s}+\sum _{k=0}^{s-1}a_{k}z^{k}\quad {\text{and}}\quad \sigma (z)=\sum _{k=0}^{s}b_{k}z^{k}.}$$ 以多项式表示的话，上述方式阶数为p的条件会变成 $${\displaystyle \rho (e^{h})-h\sigma (e^{h})=O(h^{p+1})\quad {\text{as }}h\to 0.}$$ 此方法具有一致性，若其阶数大于等于1，也就是${\displaystyle \rho (1)=0}$ and ${\displaystyle \rho '(1)=\sigma (1)}$。

稳定性和收敛

单步法数值解会和其初始条件${\displaystyle y_{0}}$有关，而s步方法的解会和其s个初始值${\displaystyle y_{0},y_{1},\ldots ,y_{s-1}}$有关。因此需考虑此数值解在其初始值有扰动的情形下，是否会稳定。针对某微分方程的线性多步法，在特定时间区间内具有零稳定性（zero-stable），若初始值大小为ε的扰动，造成该时间区间数值解的变化维持在Kε以内，其中K为常数，不会随h步阶而变化。这称为零稳定性，因为可以用微分方程${\displaystyle y'=0}$来确认此特性(Süli & Mayers 2003，第332页)。

若特征方程ρ的根，其大小都小于等于1，且大小为1的根没有重根，则称其为满足root condition。线性多步法具有零稳定性，若且唯若其满足root condition(Süli & Mayers 2003，第335页)。

现在假如一个有一致性的线性多步法，用在够光滑的微分方程，随著${\displaystyle h\to 0}$，其起始值${\displaystyle y_{1},\ldots ,y_{s-1}}$全部收敛于起始值${\displaystyle y_{0}}$。在${\displaystyle h\to 0}$时，此数值解收敛于精确解，若且唯若其方法是零稳定的。此结果称为达尔奎斯特等价定理（Dahlquist equivalence theorem），得名自杰蒙德·达尔奎斯特（英语：Germund Dahlquist），此定理在其理念上和有限差分法的拉克斯等价定理类似。若此方法为p阶方法，则全域误差（英语：global truncation error）（特定时间内，数值解和精确解的差距）会是${\displaystyle O(h^{p})}$ (Süli & Mayers 2003，第340页)。

再者，考虑一个收敛的方法，此方法为强稳定（strongly stable）若${\displaystyle z=1}$是其中唯一大小为1的根。若此方法收敛，所有大小为1的根都不是重根，但存在不止一个大小为1的根，此方法称为相对稳定（relatively stable）。若此方法要收敛，1一定要是其中的一个根。收敛的线性多步法若不是强稳定，就会是相对稳定。

要评估刚性方程线性多步法的性能，考虑线性测试方程y' = λy。用在此方程，步阶为h的多步法，会得到线性递回关系式，其特征方程为$${\displaystyle \pi (z;h\lambda )=(1-h\lambda \beta _{s})z^{s}+\sum _{k=0}^{s-1}(\alpha _{k}-h\lambda \beta _{k})z^{k}=\rho (z)-h\lambda \sigma (z).}$$ 此方程为多步法的稳定性方程（stability polynomial）。若此方程所有根的大小都小于1，此多步法的数值解会收敛到0，此多步法会称为绝对稳定，对于hλ的值。此方法称为A稳定，若针对所有实部为负的hλ，都是绝对稳定。其绝对稳定区域是所有使多步法绝对稳定的hλ的集合(Süli & Mayers 2003，第347 & 348页)。

第一和第二达尔奎斯特壁垒

这二个结果是由杰蒙德·达尔奎斯特所证明的结果，说明了线性多步法收敛阶数和A稳定性上的重要限制。第一达尔奎斯特壁垒是在Dahlquist (1956)所证明，第二壁垒则是在Dahlquist (1963)证明。

第一达尔奎斯特壁垒

第一达尔奎斯特壁垒指出：零稳定且线性的q步多步法，若q是奇数，其收敛阶数无法超过q + 1，若q是偶数，其收敛阶数无法超过q + ２。若此方法是显式方法，其收敛阶数无法超过q(Hairer, Nørsett & Wanner 1993，Thm III.3.5)。

第二达尔奎斯特壁垒

第二达尔奎斯特壁垒指出：显式线性多段法不会是A稳定。而且，A稳定隐式线性多段法的最大阶数是2。在收敛阶数为2的A稳定线性多段法中，梯形法的误差常数最小(Dahlquist 1963，Thm 2.1 and 2.2)。

相关条目

• 能量漂移（英语：Energy drift）

参考资料

• Bashforth, Francis, An Attempt to test the Theories of Capillary Action by comparing the theoretical and measured forms of drops of fluid. With an explanation of the method of integration employed in constructing the tables which give the theoretical forms of such drops, by J. C. Adams, Cambridge, 1883.
• Butcher, John C., Numerical Methods for Ordinary Differential Equations, John Wiley, 2003, ISBN 978-0-471-96758-3.
• Dahlquist, Germund, Convergence and stability in the numerical integration of ordinary differential equations, Mathematica Scandinavica, 1956, 4: 33–53, doi:10.7146/math.scand.a-10454 .
• Dahlquist, Germund, A special stability problem for linear multistep methods, BIT, 1963, 3: 27–43, ISSN 0006-3835, S2CID 120241743, doi:10.1007/BF01963532.
• Goldstine, Herman H., A History of Numerical Analysis from the 16th through the 19th Century, New York: Springer-Verlag, 1977, ISBN 978-0-387-90277-7.
• Hairer, Ernst; Nørsett, Syvert Paul; Wanner, Gerhard, Solving ordinary differential equations I: Nonstiff problems 2nd, Berlin: Springer Verlag, 1993, ISBN 978-3-540-56670-0.
• Hairer, Ernst; Wanner, Gerhard, Solving ordinary differential equations II: Stiff and differential-algebraic problems 2nd, Berlin, New York: Springer-Verlag, 1996, ISBN 978-3-540-60452-5.
• Iserles, Arieh, A First Course in the Numerical Analysis of Differential Equations, Cambridge University Press, 1996, Bibcode:1996fcna.book.....I, ISBN 978-0-521-55655-2.
• Milne, W. E., Numerical integration of ordinary differential equations, American Mathematical Monthly (Mathematical Association of America), 1926, 33 (9): 455–460, JSTOR 2299609, doi:10.2307/2299609.
• Moulton, Forest R., New methods in exterior ballistics, University of Chicago Press, 1926.
• Quarteroni, Alfio; Sacco, Riccardo; Saleri, Fausto, Matematica Numerica, Springer Verlag, 2000, ISBN 978-88-470-0077-3.
• Süli, Endre; Mayers, David, An Introduction to Numerical Analysis, Cambridge University Press, 2003, ISBN 0-521-00794-1.

外部链接

• 埃里克·韦斯坦因. Adams Method. MathWorld.