常微分方程的数值方法

线性多步法

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

线性多步法(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=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[&\frac73f(x_n,y(x_n))-\frac23f(x_{n-1},y(x_{n-1}))\\ &+\frac13f(x_{n-2},y(x_{n-2}))] \end{aligned} \]

这样，有格式

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

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

\[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 \]

或者，取$y_{n-1}=y(x_{n-1})$, $f(x_j,y_j)=f(x_j,y(x_j))=y'(x_j)$，然后做Taylor展开即可。 可以得到，误差为

\[T_{n+1}=\frac13h^4y^{(4)}(\xi) \]

格式为3阶3步的格式。

例 2. 建立$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} \]

线性多步法

通常的$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} \]

其中$f_j=f(x_j,y_j)$, $\alpha$, $\beta$为系数。

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

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

\[LTE=y(x_{i+1})-(\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}))) \]

利用$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{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} \]

待定系数法构造

命题 .当$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$

例 3. 当k=1时，可以得到条件

\[\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)$

例 4. 推导形如

\[\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} \]
\[\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个方程。有无穷组解。

线性多步法的使用

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

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

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

线性多步法的相容性

当线性多步格式

\[\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$的高阶无穷小量时，格式是相容（Consistency）的。利用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 \]

线性多步法的稳定性

通常的$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}$ 为该$k$步法的特征多项式。设$\lambda_1$, $\lambda_2$, $\cdots$, $\lambda_k$是特征方程$P(\lambda)=0$的根。 则有

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

目录

本节读完

例 5.

[#ex9-1-0].