线性多步法

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

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

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

## 定义

常微分方程的数值方法设法找到以下初值问题的近似解$$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。