差分方程的绝对稳定性

常微分方程的数值方法

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

差分方程的绝对稳定性

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

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

解. 方程有真解

\[y(x)=e^{-30x} , \quad x\in[0, 0.5] \]

是一个单调递减趋于0的函数。

取步长$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}=\frac1{1+30h}y_i=\frac14y_i, y_0=1$

注. 这里， 隐格式直接解出来，没有使用预估-校正过程

上述格式如果用预估-校正过程，会得到

\[\begin{aligned} \tilde y_{i+1}=&y_i+hf(x_i,y_i)=y_i-30hy_i=(1-30h)y_i \\ y_{i+1}=&y_i+hf(x_{i+1},\tilde y_{i+1})=y_i-30h\tilde y_{i+1} \end{aligned} \]

因而

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

改进的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}$

num-ode-stable

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

  • 在向前Euler格式中，取$h=0.01$，则有$y_{i+1}=0.7y_i$，同样可以得到稳定的解。
  • 也就是说， 向前Euler格式的稳定性与步长的大小有关

num-ode-stable-2

$h=0.05$

num-ode-stable-3

$h=0.025$

更一般地，对于如下方程

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

下面，分析一下舍入误差的作用。

考虑最简单的情形： 初值带有误差$e_0$， 以后的计算中都没有舍入误差。

假定初值$y_0$带了误差$e_0$后，变成了$z_0$，则计算过程变成了，

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

与$y_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$。这个区域称为稳定区域。

简单来说就是，是否存在步长$h$，使得格式在计算$y'=\lambda y$时，能够得到稳定的数值结果。

单步格式的绝对稳定性

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

解. 作用到微分方程$y'(x)=\lambda y$后，得到

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

因而

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

稳定区域为

\[|1+\lambda h|<1 \]

从复平面上看，$\lambda h$是一个以$(-1,0)$为圆心，1为半径的圆。

例 3. 向后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格式是绝对稳定的。

例 4. 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) \]

稳定区域满足

\[|1+\lambda h+\frac12(\lambda h)^2+\frac16(\lambda h)^3|<1 \]

Runge-Kutta法的绝对稳定区间

误差放大率 稳定区间
1阶 $1+\lambda h$ (-2,0)
2阶 $1+\lambda h+\frac1{2!}(\lambda h)^2$ (-2,0)
3阶 $1+\lambda h+\frac12(\lambda h)^2+\frac16(\lambda h)^3$ (-2.51,0)
4阶 $1+\lambda h+\frac12(\lambda h)^2+\frac16(\lambda h)^3+\frac1{24}(\lambda h)^4$ (-2.78,0)

ode-stable-rk4

多步格式的绝对稳定性

例 5. 中心差商格式的绝对稳定性

解. 中心差商格式应用到$y'=\lambda y$

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

因而，误差满足方程

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

这是一个多步的公式。考虑初值$e_0$, $e_1$对以后的影响。

方程对应的特征方程是

\[\xi^2-\lambda h \xi-1=0 \]

它的根

\[\xi_{1,2}=\lambda h \pm \sqrt{1+(\lambda h)^2} \]

因而，误差有公式

\[e_{i}=a\xi_1^i+b\xi_2^i \]

其中$a$, $b$由$e_0$, $e_1$给出，即满足

\[\begin{cases} e_0=a+b \\ e_1=a\xi_1+b\xi_2 \end{cases} \]

由$\lambda h<0$，则

\[|\xi_2|=|\lambda h-\sqrt{1+(\lambda h)^2}|>|\sqrt{1+(\lambda h)^2}|>1 \]

因而，中心差商格式对于任意$h>0$都不稳定。

中心差商计算$y'(x)=-30y$，区间$[0,0.5]$。得到不同区间数下的结果 

 5 Error:  3038.000000305902
10 Error:  34388.4999996941 
20 Error:  104857.60000055241
40 Error:  74742.79664809995

40个区间的结果 

[1.0 0.625 0.53125 0.2265625 0.361328125 -0.04443359375 0.3946533203125
 -0.340423583984375 0.6499710083007812 -0.8279018402099609
 1.270897388458252 -1.78107488155365 2.6067035496234894 -3.736102543771267
 5.40878045745194 -7.792687886860222 11.253296372597106 -16.23266016630805
 23.427791497328144 -33.80350378930416 48.78041933930626
 -70.38881829378386 101.57203305964416 -146.56784308851698
 211.4979153760319 -305.1912796205409 440.39137509143757
 -635.4848109391191 917.004983295777 -1323.2385484109518
 1909.4338946039911 -2755.313969363945 3975.91937162695 -5737.253498084157
 8278.859495190069 -11946.39811947671 17238.6580847976 -24875.391683074908
 35895.20184710379 -51796.79306840275 74742.79664840584]

例 6. 2阶Adams-Bashforth格式的绝对稳定性

解. 将格式应用到方程$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 \]

ode-stable-ab2

例 7. Simpson方法的绝对稳定性

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

\[\begin{aligned} y_{i+1}=&y_{i-1}+\frac{h}3\left(f(x_{i+1},y_{i+1})+4f(x_i,y_i)+f(x_{i-1},y_{i-1})\right) \\ =&y_{i-1}+\frac{\lambda h}3\left(y_{i+1}+4 y_i+y_{i-1}\right) \end{aligned} \]

得到

\[\left(1-\frac{\lambda h}3\right)y_{i+1}-\frac{4\lambda h}3 y_i -\left(1+\frac{\lambda h}3\right)y_{i-1}=0 \]

它的特征方程为

\[\left(1-\frac{\lambda h}3\right)x^2-\frac{4\lambda h}3 x -\left(1+\frac{\lambda h}3\right)=0 \]

特征根为

\[\xi_{1,2}=\frac{{\frac23}z\pm\sqrt{1+\frac13z^2}}{1-\frac13z}, z=\lambda h \]

注意到$z=\lambda h<0$，可以得到

\[|\xi_2|=\left|\frac{{\frac23}z-\sqrt{1+\frac13z^2}}{1-\frac13z}\right|>1 \]

因而，格式不是绝对稳定的。

注. 记$x=-z>0$，则

\[\begin{aligned} |\xi_2|=&\left|\frac{-{\frac23}x-\sqrt{1+\frac13x^2}}{1+\frac13x}\right| =&\frac{{\frac23}x+\sqrt{1+\frac13x^2}}{1+\frac13x} \end{aligned} \]

由$x>0$时，有

\[{\frac23}x+\sqrt{1+\frac13x^2}>{\frac23}x+1>{1+\frac13x} \]

从而，$|\xi_2|>1$

  • 可以看到，中心差商方法和Simpson方法都是稳定的格式，但不是绝对稳定的格式。 也就是说，两种格式不能在实际应用中使用，但可以用在如起步计算、预估步的计算中。
  • 一般来说，隐格式的绝对稳定性比显格式的要好，如：向后Euler格式与向前Euler格式比较。

谢谢

vertical slide 2

目录

本节读完

例 8.

8.