線性多步法

線性多步法是一種初值問題的常微分方程數值方法。在概念上，初值問題的數位方法是由初值開始，在時間上前進一小段，求解下一個點的數值。重覆此步驟，直到要求解的區間都完成為止。

單步法（像歐拉方法）只用前一個點以及其微分來計算目前的值，單步法的龍格-庫塔法會用前一個點以及和中間點中間的值來計算目前的值，不過在計算下一個值時，就不考慮之前的資訊。

多步法和單步法的差異是，多步法會考慮以前的資料，以提昇效率。因此。多步法會參考之前的數個點以及數點微分值。線性多步法會使用之前的點和微分值的線性組合。

定義

常微分方程的數值方法設法找到以下初值問題的近似解 $${\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.