常微分方程的数值方法

差分方程的绝对稳定性

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

差分方程的绝对稳定性

例 1. 分别用向前Euler，向后Euler，改进的Euler格式解初值问题

\[\begin{cases} y'(x)=-30y \\ y(0)=1 \end{cases} x\in[0,0.5] \]

解. 取步长$h=0.1$，则有

向前Euler

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

则$y_{i+1}=-2y_i, y_0=1$

向后Euler

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

则$y_{i+1}=\frac14y_i, y_0=1$

改进的Euler

\[y_{i+1}=\frac12(y_i+(1-30h)^2y_i) \]

则$y_{i+1}=\frac52y_i, y_0=1$

$x_i$ 向前Euler 向后Euler 改进的Euler 真解
0.1 -2 1/4 2.5 $4.9797\times 10^{-2}$
0.2 4 $(1/4)^2$ $2.5^2$ $2.4788\times 10^{-3}$
0.3 -8 $(1/4)^3$ $2.5^3$ $1.2341\times 10^{-4}$
0.4 16 $(1/4)^4$ $2.5^4$ $6.1442\times 10^{-6}$
0.5 -32 $(1/4)^5$ $2.5^5$ $3.0590\times 10^{-7}$

问题. 向前Euler格式、改进的Euler格式得到的结果并不稳定，为什么？

在向前Euler格式中，取$h=0.01$，则有$y_{i+1}=0.7y_i$，同样可以得到稳定的解。 事实上，取$h$足够小（$h<\frac1{30}$）后，向前Euler格式是稳定的。也就是说， 向前Euler格式的稳定性与步长的大小有关

对于如下方程

\[\begin{cases} y'(x)=\lambda y (Re(\lambda<0) \\ y(a)=y_0 \end{cases} \]

应用向前Euler格式，有

\[y_{i+1}=y_i+\lambda h y_i=(1+\lambda h)y_i \]

假定初值$y_0$带了误差后，变成了$z_0$，则有

\[z_{i+1}=(1+\lambda h)z_i \]

则有

\[y_{i+1}-z_{i+1}=(1+\lambda h)(y_i-z_i) \]

即

\[e_{i+1}=(1+\lambda h)e_i=(1+\lambda h)^{i+1}e_0 \]

可以看到，当$|1+\lambda h|<1$时，误差是递减的。也就是说，步长足够小后，向前Euler格式是稳定的。

定义 1.
差分方法称为绝对稳定(A-Stable)的，若差分方法作用到微分方程

\[y'(x)=\lambda y, Re(\lambda)<0 \]

时，对任意的初值，总存在左半复平面上的一个区域，当$\lambda h$在这个区域时，差分方程的解趋于$0$。这个区域称为稳定区域。

注. 对向前Euler格式来说，稳定区域是以$-1$为中心，$1$为半径的圆。

例 2. 向后Euler格式的绝对稳定性

解. 将向后Euler格式应用到$y'(x)=\lambda y$后，有

\[y_{i+1}=y_i+\lambda h y_{i+1} \]

所以有

\[(1-\lambda h)e_{i+1}=e_i \]

即当$\frac{1}{|1-\lambda h|}<1$为稳定区域。注意到$Re(\lambda)<0$，整个左半复平面都是稳定区域。 即，对任意$h>0$，向后Euler格式是绝对稳定的。

例 3. 3阶Runge-Kutta格式的绝对稳定性

解. 将格式应用到方程$y'(x)=\lambda y$后，有

\[\begin{cases} y_{i+1}&=y_i+\frac{h}6(k_1+4k_2+k_3) \\ k_1&=\lambda y_i \\ k_2&=\lambda y_i(1+\frac12\lambda h) \\ k_3&=\lambda y_i(1+\lambda h+(\lambda h)^2) \end{cases} \]

即有

\[y_{i+1}=y_i(1+\lambda h+\frac12(\lambda h)^2+\frac16(\lambda h)^3) \]

例 4. 2阶Adams-Bashforth格式的A稳定性

解. 将格式应用到方程$y'(x)=\lambda y$后，有

\[y_{i+1}=y_i+h(\frac32\lambda y_i-\frac12\lambda y_{i-1}) =(1+\frac32\lambda h)y_i-\frac12\lambda hy_{i-1} \]

格式要稳定，需要上式的特征多项式

\[\phi(x)=x^2-(1+\frac32\lambda h)x+\frac12\lambda h \]

的所有根的模小于1。即有

\[\left| \tfrac12 \Big(1 + \tfrac32 z \pm\sqrt{1 + z + \tfrac94 z^2}\Big) \right| < 1 , z=\lambda h \]

谢谢

vertical slide 2

目录

本节读完

例 5.

[#ex9-1-0].