線性多步法

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

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

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

## 定義

常微分方程的數值方法設法找到以下初值問題的近似解$$y' = f(t,y), \quad y(t_0) = y_0.$$

結果是在各離散時間$t_i$下，$y(t)$的近似值：

$$y_i \approx y(t_i) \quad\text{where}\quad t_i = t_0 + i h,$$ 其中$h$是時間步長（有時會表示為$\Delta t$），而$i$是整數。

多步法會用以前$s$步以內的資訊來計算下一個值。而線性多步法會用$y_i$和$f(t_i,y_i)$的線性組合，來計算$y$。因此，線性多步法可以表示如下$$\begin{align} & 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_jy_{n+j} = h\sum_{j=0}^sb_jf(t_{n+j},y_{n+j}), \end{align}$$ 而$a_s=1$。係數$a_0, \dotsc, a_{s-1}$和$b_0, \dotsc, b_s$決定所使用的方法。設計者一方面要找到精確解的理想近似解，另一方面也要設計方便計算的方法。為了簡化計算，中許多係數會是0。

可以用係數來區分顯式和隱式方法。若$b_s = 0$，此方法則為顯式法，因為可以直接計算$y_{n+s}$。若$b_s \ne 0$，則$y_{n+s}$的值和$f(t_{n+s}, y_{n+s})$有關，需設法求解方程才能得到$y_{n+s}$, 此方法為隱式法。求解方程時，常會使用迭代法（例如牛頓法）求解。

有時，顯式的多步法會用來「預估」$y_{n+s}$的值，之後再用隱式公式中來得到「校正」後的值，這稱為預估-校正方法。

## 例子

考慮以下問題$$y' = f(t,y)=y, \quad y(0) = 1.$$ 其精確解是$y(t) = e^t$。

### 單步歐拉法

歐拉法是一種簡易的單步法：

$$y_{n+1} = y_n + hf(t_n, y_n).$$ 可以將歐拉法視為是退化的顯式多步法，其步數只有一步。

用此方法，配合步長$h = \tfrac{1}{2}$求解問題$y' = y$，得到以下的解：

$$\begin{align} <pre> 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. </pre> \end{align}$$

### 二步Adams–Bashforth法

二步的Adams–Bashforth法是簡單的多步法$$y_{n+2} = y_{n+1} + \tfrac{3}{2} hf(t_{n+1},y_{n+1}) - \tfrac{1}{2} hf(t_n,y_n).$$ 此方法會用到二個值，$y_{n+1}$和$y_n$來計算下一個值$y_{n+2}$。不過，初始問題只會有一個值$y_0 = 1$。一種可能作法是用歐拉法計算$y_1$，當作第二個值。在此選擇下，Adams–Bashforth法得到結果如下（四捨五入到小數第四位）：

$$\begin{align} <pre> 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. </pre> \end{align}$$ 在$t = t_4 = 2$的解是$e^2 = 7.3891\ldots$，二步的Adams–Bashforth法比歐拉法準確。一般而言，只要步長夠小，二步的Adams–Bashforth法都會比歐拉法準確。

## 多步法分類

常見的多步法可以分為三類：Adams–Bashforth法、Adams–Moulton法及反向微分公式（BDF）。

### Adams–Bashforth法

Adams–Bashforth法是顯式法，係數是$a_{s-1}=-1$，$a_{s-2} = \cdots = a_0 = 0$，而選擇$b_j$使此方法可以到s階（因此在同一階下，此方法的係數唯一）。

s = 1, 2, 3, 4, 5的Adams–Bashforth法如下，其中的第一式就是前向歐拉法：：

$$\begin{align} 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{align}$$

係數可以用以下方法計算。用多項式插值找到$s-1$階的多項式p使得$$p(t_{n+i}) = f(t_{n+i}, y_{n+i}), \qquad \text{for } i=0,\ldots,s-1.$$ 多項式內插的拉格朗日插值法可以得到$$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\ne j}^{s-1} (t-t_{n+i}).$$ 多項式p是要求解的微分方程$y' = f(t,y)$等號右邊的局部好近似值，因此改為考慮方程$y' = p(t)$，此方法有精確解，就是p的積分。因此得到下式$$y_{n+s} = y_{n+s-1} + \int_{t_{n+s-1}}^{t_{n+s}} p(t)\,\mathrm dt.$$ 當上述公式的p取代之後，就是Adams–Bashforth法。其係數$b_j$可以由下式求得$$b_{s-j-1} = \frac{(-1)^j}{j!(s-j-1)!} \int_0^1 \prod_{i=0 \atop i\ne j}^{s-1} (u+i) \,\mathrm du, \qquad \text{for } j=0,\ldots,s-1.$$ 將$f(t, y)$用內插的p取代，引進了數量級為hs的誤差，而s步的Adams–Bashforth法確實是s階Iserles。

Adams–Bashforth法是由約翰·柯西·亞當斯求解毛細現象建模的微分方程而來，此微分方程出自弗朗西斯·巴什福斯。發表了他的理論以及亞當斯的數值方法Goldstine。

### Adams–Moulton法

Adams–Moulton法和Adams–Bashforth法類似，也是$a_{s-1} = -1$，$a_{s-2} = \cdots = a_0 = 0$，也是選擇係數b到最高可能階數。但Adams–Moulton法是隱式法。去除了$b_s = 0$的限制，s步的Adams–Moulton法可以到$s+1$階，而s步的Adams–Bashforth法只有s階。

s = 0, 1, 2, 3, 4的Adams–Moulton法如下，其中前兩個分別是反向歐拉法和梯形法則： respectively:

$$\begin{align} y_{n} &= y_{n-1} + h f(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{align}$$

Adams–Moulton法的推導類似Adams–Bashforth法，不過，用的內插法不只用$t_{n-1},\dots, t_{n-s}$的點，也會用$t_n$點，其係數為$$b_{s-j} = \frac{(-1)^j}{j!(s-j)!} \int_0^1 \prod_{i=0 \atop i\ne j}^{s} (u+i-1) \,\mathrm du, \qquad \text{for } j=0,\ldots,s.$$

Adams–Moulton法是由約翰·柯西·亞當斯獨立完成，類似Adams–Bashforth法。福里斯特·雷·莫爾頓的名字出現在方法，是因為他發現此方法可以和Adams–Bashforth法一起使用，成為一對預估-校正方法Moulton。也有類似的想法，亞當斯是用牛頓法求解隱式方程Hairer。

### 反向微分公式（BDF）

反向微分公式是$b_{s-1} = \cdots = b_0 = 0$的隱式法，選擇其他係數使階數可以到s階（最大可能階數）。此方法特別用來求得剛性方程。

## 分析

微分方程數值方法分析的中心概念，包括有收斂、階數和穩定性。

### 收斂、階數

第一個問題是此方法是否具有一致性：以下差分方程$$\begin{align} & 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{align}$$ 是否是$y' = f(t,y)$微分方程的良好近似？更準確的說法，多步法具有一致性，當其步階h趨近於零時，局部截尾誤差趨近於零的速度比h還快，其中局部截尾誤差定義為此方法的結果$y_{n+s}$，和時間$t_{n+s}$時方程精確解的誤差（假設所有之前的值$y_{n+s-1}, \ldots, y_n$都是正確的）。利用泰勒級數的計算可得線性多步法具有一致性，若且唯若$$\sum_{k=0}^{s-1} a_k = -1 \quad\text{and}\quad \sum_{k=0}^s b_k = s + \sum_{k=0}^{s-1} k a_k.$$

上述所有的方式都具有一致性Hairer。

若方法具有一致性，下一個問題是：近似微分方程的差分方程，可以近似到多好？多步法其階數為p若在其步階h趨近於零時，其局部誤差約為$O(h^{p+1})$。這和以下方法係數的條件等效：

$$\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方法的階數為$s+1$Hairer。

這些條件可以用特徵多項式表示：

$$\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的條件會變成$$\rho(e^h) - h\sigma(e^h) = O(h^{p+1}) \quad \text{as } h\to 0.$$ 此方法具有一致性，若其階數大於等於1，也就是$\rho(1)=0$ and $\rho'(1)=\sigma(1)$。

### 穩定性和收斂

單步法數值解會和其初始條件$y_0$有關，而s步方法的解會和其s個初始值$y_0, y_1, \ldots, y_{s-1}$有關。因此需考慮此數值解在其初始值有擾動的情形下，是否會穩定。針對某微分方程的線性多步法，在特定時間區間內具有零穩定性（zero-stable），若初始值大小為ε的擾動，造成該時間區間數值解的變化維持在Kε以內，其中K為常數，不會隨h步階而變化。這稱為零穩定性，因為可以用微分方程$y' = 0$來確認此特性Süli。

若特徵方程ρ的根，其大小都小於等於1，且大小為1的根沒有重根，則稱其為滿足root condition。線性多步法具有零穩定性，若且唯若其滿足root conditionSüli。

現在假如一個有一致性的線性多步法，用在夠光滑的微分方程，隨著$h \to 0$，其起始值$y_1, \ldots, y_{s-1}$全部收斂於起始值$y_0$。在$h \to 0$時，此數值解收斂於精確解，若且唯若其方法是零穩定的。此結果稱為達爾奎斯特等價定理（Dahlquist equivalence theorem），得名自傑蒙德·達爾奎斯特，此定理在其理念上和有限差分法的拉克斯等價定理類似。若此方法為p階方法，則全域截斷誤差（特定時間內，數值解和精確解的差距）會是$O(h^p)$ Süli。

再者，考慮一個收斂的方法，此方法為強穩定（strongly stable）若$z=1$是其中唯一大小為1的根。若此方法收斂，所有大小為1的根都不是重根，但存在不止一個大小為1的根，此方法稱為相對穩定（relatively stable）。若此方法要收斂，1一定要是其中的一個根。收斂的線性多步法若不是強穩定，就會是相對穩定。

要評估剛性方程線性多步法的性能，考慮線性測試方程y = λy。用在此方程，步階為h的多步法，會得到線性遞迴關係式，其特徵方程為$$\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。

## 第一和第二達爾奎斯特壁壘

這二個結果是由傑蒙德·達爾奎斯特所證明的結果，說明了線性多步法收斂階數和A穩定性上的重要限制。第一達爾奎斯特壁壘是在所證明，第二壁壘則是在證明。

### 第一達爾奎斯特壁壘

第一達爾奎斯特壁壘指出：零穩定且線性的q步多步法，若q是奇數，其收斂階數無法超過q + 1，若q是偶數，其收斂階數無法超過q + ２。若此方法是顯式方法，其收斂階數無法超過qHairer。

### 第二達爾奎斯特壁壘

第二達爾奎斯特壁壘指出：顯式線性多段法不會是A穩定。而且，A穩定隱式線性多段法的最大階數是2。在收斂階數為2的A穩定線性多段法中，梯形法的誤差常數最小Dahlquist。