线性多步法

常微分方程的数值方法

张瑞
中国科学技术大学数学科学学院

线性多步法(Linear multistep method)

基于数值积分的构造方法

对微分方程在区间$[x_{n-p},x_{n+1}]$上积分，有

\[\int_{x_{n-p}}^{x_{n+1}}y'(x)dx=\int_{x_{n-p}}^{x_{n+1}}f(x,y(x))dx \]

即有

\[y({x_{n+1}})-y({x_{n-p}})=\int_{x_{n-p}}^{x_{n+1}}f(x,y(x))dx \]

这样，用数值积分来近似$\displaystyle\int_{x_{n-p}}^{x_{n+1}}f(x,y(x))dx$即可得到不同的格式。数值积分的误差即为格式的局部截断误差。

  • 以$x_{n-q},x_{n-q+1},\cdots,x_n$为积分节点得到的格式为显格式，

    \[\begin{aligned} \int_{x_{n-p}}^{x_{n+1}}f(x,y(x))=&h\beta_0f(x_n,y(x_n))+h\beta_1f(x_{n-1},y(x_{n-1})) \\ &+\cdots+h\beta_qf(x_{n-q},y(x_{n-q})) \\ &+hT \end{aligned} \]

    其中$h$是步长，$h\beta_i$是积分系数，$hT$为局部截断误差。 积分节点共有$q+1$个，格式是至少$q+1$阶的多步格式。

  • 以$x_{n-q+1},x_{n-q+2},\cdots,x_{n+1}$为积分节点得到的格式为隐格式。

    \[\begin{aligned} \int_{x_{n-p}}^{x_{n+1}}f(x,y(x))\approx &h\beta_{-1}f(x_{n+1},y(x_{n+1}))+h\beta_0f(x_{n},y(x_{n})) \\ &+\cdots+h\beta_{q-1}f(x_{n-q+1},y(x_{n-q+1})) \\ \end{aligned} \]

例 1. 在等距步长为$h$的分割中，建立$p=0$, $q=1$的隐格式。

解. $p=0$，表示积分区间是$\displaystyle\int_{x_n}^{x_{n+1}}$，

\[y(x_{n+1})=y(x_n)+\int_{x_n}^{x_{n+1}}f(x,y)dx \]

$q=1$的隐格式，表示积分节点是$x_{n+1}$, $x_n$共2个。 这是梯形积分公式

\[\begin{aligned} y(x_{n+1})=y(x_n)+\frac{h}2(f(x_{n+1}, y(x_{n+1}))+f(x_n, y(x_n))) \\ -\frac{f''(\xi_n)}{12}h^3 \end{aligned} \]

得到格式为

\[y_{n+1}=y_n+\frac{h}2(f(x_n,y_n)+f(x_{n+1},y_{n+1})) \]

称为梯形格式。

局部截断误差是$-\frac{f''(\xi_n)}{12}h^3=O(h^3)$。

因此梯形公式是2阶的单步隐格式。

例 2. 在等距步长为$h$的分割中，建立$p=1$, $q=2$的显格式

解. $p=1$，则积分区间为$\displaystyle\int_{x_{n-1}}^{x_{n+1}}$，

$q=2$的显示格式，则积分节点为$x_n,x_{n-1},x_{n-2}$

\[\begin{aligned} h\beta_0=&\int_{x_{n-1}}^{x_{n+1}}\frac{(x-x_{n-1})(x-x_{n-2})}{(x_n-x_{n-1})(x_n-x_{n-2})}dx=\frac73h \\ h\beta_1=&\int_{x_{n-1}}^{x_{n+1}}\frac{(x-x_{n})(x-x_{n-2})}{(x_{n-1}-x_{n})(x_{n-1}-x_{n-2})}dx=-\frac23h \\ h\beta_2=&\int_{x_{n-1}}^{x_{n+1}}\frac{(x-x_{n})(x-x_{n-1})}{(x_{n-2}-x_{n})(x_{n-2}-x_{n-1})}dx=\frac13h \\ \end{aligned} \]

可以得到

\[\begin{aligned} y(x_{n+1})-y(x_{n-1})\approx h\bigg[&\frac73f(x_n,y(x_n))-\frac23f(x_{n-1},y(x_{n-1}))\\ &+\frac13f(x_{n-2},y(x_{n-2}))\bigg] \end{aligned} \]

这样，有格式

\[y_{n+1}=y_{n-1}+h\left(\frac73f(x_n,y_n)-\frac23f(x_{n-1},y_{n-1})+\frac13f(x_{n-2},y_{n-2})\right) \]

局部截断误差： 为数值积分的误差

\[hT=\int_{x_{n-1}}^{x_{n+1}}\frac{y^{(4)}(\xi)}{3!}(x-x_n)(x-x_{n-1})(x-x_{n-2})dx \]

用Taylor展开的方法算局部截断误差

当$y_n$, $y_{n-1}$, $y_{n-2}$是精确时，有

\[\begin{aligned} y_{n-1}=& y(x_{n-1}),\\ f(x_j,y_j)=&f(x_j,y(x_j))=y'(x_j), j=n,n-1,n-2 \end{aligned} \]

代入格式后，得到

\[y_{n+1}=y(x_{n-1})+\frac{h}3(7y'(x_n)-2y'(x_{n-1})+y'(x_{n-2})) \]

上式在$x_{n-1}$处做Taylor展开，

\[\begin{aligned} y'(x_n)=y'(x_{n-1})+hy''(x_{n-1})+\frac{y'''(x_{n-1})}{2!}h^2 \\ +\frac{y^{(4)}(x_{n-1})}{3!}h^3+o(h^3) \end{aligned} \]
\[\begin{aligned} y'(x_{n-2})=y'(x_{n-1})-hy''(x_{n-1})+\frac{y'''(x_{n-1})}{2!}h^2 \\ -\frac{y^{(4)}(x_{n-1})}{3!}h^3+o(h^3) \end{aligned} \]

因而，

\[\begin{aligned} y_{n+1}=&y(x_{n-1})+2hy'(x_{n-1})+\frac{4h^2}{2}y''(x_{n-1}) \\ &+\frac{8h^3}{6}y'''(x_{n-1})+\frac{h^4}{3}y^{(4)}(x_{n-1})+o(h^4) \end{aligned} \]

可以得到，误差为

\[\begin{aligned} T_{n+1}=&\left(\frac{(2h)^4}{4!}-\frac13h^4\right)y^{(4)}(x_{n-1})+o(h^4)\\ =&\frac13h^4y^{(4)}(x_{n-1})= O(h^4) \end{aligned} \]

格式为3阶3步的格式。

例 3. 建立$p=2$, $q=2$的隐格式

解. $p=2$,则积分区间为$\displaystyle\int_{x_{i-2}}^{x_{i+1}}$

$q=2$的隐格式，则积分节点为$x_{i+1}$,$x_i$, $x_{i-1}$

可以得到，积分系数为

\[\begin{aligned} h\beta_{-1}=&\frac34h \\ h\beta_0=&0 \\ h\beta_{1}=&\frac94h \end{aligned} \]

格式为

\[y_{i+1}=y_{i-2}+h(\frac34f(x_{n+1},y_{n+1})+\frac94f(x_{n-1},y_{n-1})) \]

局部截断误差

\[T_{i+1}=-\frac38h^4y^{(4)}(\xi) \]

格式是3步3阶隐格式

定义 1.
$p=0$时的公式，称为Adams公式

以下，简记$f_j=f(x_j,y_j)$，

  • 二阶显式Adams公式
    \[y_{n+1}=y_n+\frac{h}2(3f(x_n,y_n)-f(x_{n-1},y_{n-1})) \]
  • 三阶显式Adams公式
    \[y_{n+1}=y_n+\frac{h}{12}(23f_n-16f_{n-1}+5f_{n-2}) \]
  • 三阶隐式Adams公式
    \[y_{n+1}=y_n+\frac{h}{12}(5f_{n+1}+8f_n-f_{n-1}) \]

利用数值积分的误差表达式，可以得到

  • Adams显式格式（Adams-Bashforth）的误差为

    \[\begin{aligned} hT=&\int_{x_n}^{x_{n+1}}\frac{y^{(q+2)}(\xi)}{(q+1)!}(x-x_n)\cdots(x-x_{n-q})dx \\ =&\frac{y^{(q+2)}(\xi)}{(q+1)!}\int_0^1t(t+1)\cdots(t+q)dt \end{aligned} \]

  • Adams隐式格式（Adams-Moulton）的误差为

    \[\begin{aligned} hT=&\int_{x_n}^{x_{n+1}}\frac{y^{(q+2)}(\xi)}{(q+1)!}(x-x_{n+1})\cdots(x-x_{n-q+1})dx \\ =&\frac{y^{(q+2)}(\xi)}{(q+1)!}\int_0^1(t-1)t\cdots(t+q-1)dt \end{aligned} \]

线性多步法

格式

\[\begin{aligned} y_{i+1}=&\alpha_0y_i+\alpha_1y_{i-1}+\cdots+\alpha_{k-1}y_{i-k+1} \\ &+h(\beta_{-1}f_{i+1}+\beta_0f_{i}+\cdots+\beta_{k-1}f_{i-k+1}) \end{aligned} \]

被称为$k$步线性多步格式，其中$f_j=f(x_j,y_j)$, $\alpha$, $\beta$为系数。

当$\beta_{-1}\neq 0$时，格式是显式的，否则是隐式的。

注. 显然，前面用数值积分方法得到的格式和Euler公式都是线性多步格式

它的局部截断误差为（假定$y_{i},\cdots,y_{i-k+1}$是精确的），

\[LTE=y(x_{i+1})-\left(\sum_{j=0}^{k-1}\alpha_jy(x_{i-j})+h\sum_{j=-1}^{k-1}\beta_jf(x_{i-j},y(x_{i-j}))\right) \]

利用$f(x_{i-j},y(x_{i-j}))=y'(x_{i-j})$，把上式在$x_i$处Taylor展开，可以得到误差表达式。

\[\begin{aligned} LTE=&y(x_i)+hy'(x_i)+\frac{h^2}{2!}y''(x_i)+\cdots \\ & -\sum_{j=0}^{k-1}\alpha_j[y(x_i)+(-jh)y'(x_i)+\frac{(-jh)^2}{2!}y''(x_i)+\cdots ]\\ & -h\sum_{j=-1}^{k-1}\beta_j[y'(x_i)+(-jh)y''(x_i)+\frac{(-jh)^2}{2!}y'''(x_i)+\cdots ] \end{aligned} \]
\[\begin{aligned} LTE=&y(x_i)+hy'(x_i)+\frac{h^2}{2!}y''(x_i)+\cdots \\ & -\sum_{j=0}^{k-1}\alpha_j[y(x_i)+(-jh)y'(x_i)+\frac{(-jh)^2}{2!}y''(x_i)+\cdots ]\\ & -h\sum_{j=-1}^{k-1}\beta_j[y'(x_i)+(-jh)y''(x_i)+\frac{(-jh)^2}{2!}y'''(x_i)+\cdots ] \end{aligned} \]

记

\[\begin{cases} c_0=&\displaystyle 1-\sum_{j=0}^{k-1}\alpha_j \\ c_m=&\displaystyle\frac1{m!}-\frac1{m!}\sum_{j=0}^{k-1}(-j)^m\alpha_j-\frac1{(m-1)!}\sum_{j=-1}^{k-1}(-j)^{m-1}\beta_j , m=1,\cdots \end{cases} \]

则有

\[\begin{aligned} LTE=&\sum_{m=0} c_m h^m y^{(m)}(x_i) \\ =&c_0 y(x_i)+c_1 h y'(x_i)+c_2 h^2 y''(x_i)+\cdots \end{aligned} \]

命题. 当$c_i=0, i=0,1,\cdots,p$, $c_{p+1}\neq 0$时，格式的局部截断误差是$O(h^{p+1})$，是$p$阶格式。

待定系数法构造

要构造$p$阶的格式，取合适的$\alpha_i$, $\beta_i$，使得

\[c_i=0, i=0,1,\cdots,p \]

例 4. 当k=1时，格式为

\[y_{i+1}=\alpha_0 y_i +h (\beta_{-1}f_{i+1}+\beta_0 y_i) \]

可以得到条件

\[\alpha_0=1, \beta_{-1}+\beta_0=1 \]
  • $\beta_{-1}=0$时，可以得到向前Euler方法 $y_{j+1}=y_j+hf_i$
  • $\beta_{-1}=1$时，可以得到向后Euler方法 $y_{j+1}=y_j+hf_{i+1}$
  • $\beta_{-1}=\beta_0=\frac12$时，可以得到梯形公式 $y_{j+1}=y_j+\frac{h}2(f_{i+1}+f_j)$

例 5. 推导形如

\[\begin{aligned} y_{i+1}=&\alpha_0y_i+\alpha_1y_{i-1}+\alpha_2y_{i-2} \\ &+h(\beta_{-1}f_{i+1}+\beta_0f_i+\beta_1f_{i-1}+\beta_2f_{i-2}) \end{aligned} \]

的隐式线性三步格式，确定系数，使格式为4阶方法。

解. 列出方程为

\[\begin{array}{l} \alpha_0+\alpha_1+\alpha_2=1 \\ \alpha_1+2\alpha_2-\beta_{-1}-\beta_0-\beta_1-\beta_2=-1 \\ \frac12\alpha_1+2\alpha_2+\beta_{-1}-\beta_1-2\beta_2=\frac12 \\ \frac16\alpha_1+\frac43\alpha_2-\frac12\beta_{-1}-\frac12-\beta_1-2\beta_2=-\frac16 \\ \frac1{24}\alpha_1+\frac23\alpha_2+\frac16\beta_{-1}-\frac16\beta_1-\frac43\beta_2=\frac1{24} \\ \frac1{120}\alpha_1+\frac4{15}\alpha_2-\frac1{24}\beta_{-1}-\frac1{24}\beta_1-\frac23\beta_2\neq-\frac1{120} \end{array} \]

7个未知数，5个方程。有无穷组解。

Simpson格式

\[y_{i+1}=y_{i-1}+\frac{h}3(f_{i+1}+4f_i+f_{i-1}) \]

Hamming格式

\[y_{i+1}=\frac18(8y_i-y_{i-2})+\frac{3h}8(f_{i+1}+2f_i-f_{i-1}) \]

4阶Adams隐格式

\[y_{i+1}=y_i+\frac{h}{24}(9f_{i+1}+19f_i-5f_{i-1}+f_{i-2}) \]
\[\begin{aligned} &\alpha_0+\alpha_1+\alpha_2=1 \\ &\alpha_1+2\alpha_2-\beta_{-1}-\beta_0-\beta_1-\beta_2=-1 \\ &\frac12\alpha_1+2\alpha_2+\beta_{-1}-\beta_1-2\beta_2=\frac12 \\ &\frac16\alpha_1+\frac43\alpha_2-\frac12\beta_{-1}-\frac12-\beta_1-2\beta_2=-\frac16 \\ &\frac1{24}\alpha_1+\frac23\alpha_2+\frac16\beta_{-1}-\frac16\beta_1-\frac43\beta_2=\frac1{24} \\ &\frac1{120}\alpha_1+\frac4{15}\alpha_2-\frac1{24}\beta_{-1}-\frac1{24}\beta_1-\frac23\beta_2\neq-\frac1{120} \end{aligned} \]

线性多步法的使用

  1. 如何起步？ 只能使用最多低一阶的格式作为起步计算
  2. 隐格式如何计算？ 只能使用最多低一阶的格式做预估步的计算

预估-校正时，若预估步的局部截断误差为$T^*_{n+1}$，校正步的局部截断误差为$T_{n+1}$， 则整个预估-校正格式的局部截断误差为

\[T_{n+1}+hT^*_{n+1}\beta_0\frac{\partial f}{\partial y}(\xi) \]

线性多步法的相容性

定义 2.
格式的相容性(Consistency)是指，当$h\to0$时，差分方程收敛到微分方程。

当线性多步格式

\[\begin{aligned} y_{i+1}=&\alpha_0y_i+\alpha_1y_{i-1}+\cdots+\alpha_{k-1}y_{i-k+1} \\ &+h(\beta_{-1}f_{i+1}+\beta_0f_{i}+\cdots+\beta_{k-1}f_{i-k+1} \end{aligned} \]

的局部截断误差是步长$h$的高阶无穷小量时，格式是相容的。

利用Taylor展开，可以证明

定理 1.
线性多步格式是相容的充要条件是

\[\sum_{j=0}^{k-1}\alpha_{j}=1 , \sum_{j=-1}^{k-1}\beta_j=k-\sum_{j=0}^{k-1}(k-1-j)\alpha_j \]

线性多步法的稳定性

定义 3.
格式的稳定性是指差分方程对舍入误差具有抑制作用（不会被无限地放大）。

一个$k$步线性多步格式

\[\begin{aligned} y_{i+1}=&\alpha_0y_i+\alpha_1y_{i-1}+\cdots+\alpha_{k-1}y_{i-k+1} \\ &+h(\beta_{-1}f_{i+1}+\beta_0f_{i}+\cdots+\beta_{k-1}f_{i-k+1} \end{aligned} \]

的特征多项式是

\[P(\lambda)=\lambda^k-\alpha_0\lambda^{k-1}-\alpha_1\lambda^{k-2}\cdots-\alpha_{k-1} \]

设$\lambda_1$, $\lambda_2$, $\cdots$, $\lambda_k$是特征方程$P(\lambda)=0$的根。 则有

定理 2.
$k$步法格式稳定的充要条件是，$|\lambda_i|\leq 1, \forall i=1,2,\cdots,k$，且模为$1$的都是单根。

例 6. 向前Euler格式

\[y_{i+1}=y_i+hf(x_i, y_i) \]

的特征多项式是

\[P(\lambda)=\lambda-1 \]

具有单根$1$，因而是稳定的。

例 7. 向后Euler格式

\[y_{i+1}=y_i+hf(x_{i+1}, y_{i+1}) \]

的特征多项式是

\[P(\lambda)=\lambda-1 \]

具有单根$1$，因而是稳定的。

例 8. 中心差商格式

\[y_{i+1}=y_{i-1}+2hf(x_i, y_i) \]

的特征多项式是

\[P(\lambda)=\lambda^2-1 \]

具有单根$1$, $-1$，因而是稳定的。

例 9. Simpson格式

\[y_{i+1}=y_{i-1}+\frac{h}3(f_{i+1}+4f_i+f_{i-1}) \]

的特征多项式是

\[P(\lambda)=\lambda^2-1 \]

具有单根$1$, $-1$，因而是稳定的。

目录

本节读完

例 10.

10.