第5节绝对稳定性和绝对稳定区域

数值稳定性问题不但与计算公式有关，而且与微分方程本身的 性质有关。 因此，需要研究当步长h固定时用方法（3.2）计算， 若初始近似或计算过程有误差,对以后计算结果的影响是否增长， 如果误差不增长这个方法就是计算稳定的，实际对于一个相容且 渐进稳定的方法，当步长h取不同时计算结果差别很大。 例5.1 用3阶显式Adams公式
它的根 λ j < 1 , j=1,2,3。 这是因为，由 94 1 3 2 max λ 2 , 1 λ = 73λ + 16λ − 5 ≤ 96 96
(
)
即可知
λ j < 1 , j=1,2,3。 故对h=0.0012, 是绝对稳定的。
而当h=0.125时，h = −
3
2204 2 1600 500 λ + λ − λ+ =0 96 96 96
∑ (α
j =0
k
j
− h β j ) en + j = ηn
（1.5.6）
如前引进k维列向量
En = ( en + k −1 , en + k − 2 , L, en ) 和 ηn = (ak−1η n , 0, L×k阶矩阵
−1 −1 −1 ⎛ −ak ak −1 − ak ak − 2 L −ak−1a1 −ak a0 ⎞ ⎜ ⎟ L 0 0 0 ⎟ ⎜ 1 (1.5.7) ⎜ ⎟ O O M M C= ⎜ ⎟ O O M ⎟ ⎜ ⎜ ⎟ 1 0 ⎝ ⎠ 其中 a j = α j − h β j , j = k , k − 1, L , 0 。假定h充分小使得ak≠0
又 u = e 满足u ′ = μ u ，若方法为p阶，则
μu
ρ (e ) − h σ (e ) =
μu
μu
∑αe
j =0 j
k
μ ( t + jh)
− h∑ βj e
j =0
k
μ ( t + jh)
p +1 O h = ( ) = O ( h p+1 )
h h p +1 ρ ( e ) − h σ ( e ) = O ( h ) 将左端分解因式为： 从而
进一步
u′ = Tv ′ = ATv , v ′ = T −1 ATv = diag ( λ1 , L , λm ) v
T −1 AT = diag ( λ1 , L , λm
)
即
vi′ = λi vi , i = 1, L , m
因此作为模型,我们考虑初值问题：
u′ = μ u
的绝对稳定性, 其中 μ 为实或复常数。
un + 3
求初值问题 u ′ = −100 u − t 3 + 3t 2 , u (0) = 0 的数值解。 步长分别取 h = 0.00125 及 h = 0.125 解 此问题的精确解为u(t)=t3 ,公式的第一特征多项式为：
h = un + 2 + (23 f n + 2 − 16 f n +1 + 5 f n ) 12
注意，
3
100 。 特征方程为: 8
1600 500 2204 2 (− 10) + × (− 10 ) − × (− 10 ) + >0 96 96 96 2204 1600 500 2 × ( −∞ ) − × ( −∞ ) + <0 ( −∞ ) + 96 96 96
3
故在(-∞,-10)中有一个根，故方法对h=0.125是不绝对稳定的。 结论：存在着（零）稳定而不是绝对稳定的情况
(
)
ρ (λ ) = λ3 − λ2 = λ2 (λ − 1)
满足根条件,是稳定的。 但当h取不同时计算结果差别很大,见表
数值结果 tn 0.000 0.125 0.250 0.375 0.500 0.625 0.750 0.875 1.000 h=0.125 un u(tn)-un 0.00000000 0.00195307 0.00000004 0.01562409 - 0.00000091 0.05275568 0.00002130 0.12449563 -0.00050437 0.25408001 0.00993938 0.18516620 -0.23670880 6.27264466 5.60272278 -131.625987 -132.625987 h=0.00125 un u(tn)-un 0.00000000 0.00195316 0.00000003 0.01562502 0.00000002 0.05273460 - 0.00000002 0.12499999 - 0.00000001 0.24414061 - 0.00000002 0.42187498 - 0.00000002 0.66992190 0.00000002 0.99999999 - 0.00000001
j j =0
k
σ (λ ) = ∑ β j λ j
j =0
k
（1.5.10）
则 C的特征方程为
ρ (λ ) − h σ (λ ) = 0
（1.5.11）
由引理1.1(ⅲ),可知(1.5.9)成立的充要条件是(1.5.11)的根都在 单位圆内。 事实上,当矩阵
ρ ( C ) < 1, 则存在一种矩阵的算子范数‖·‖,使得 C < 1。 则由
对一般非线性常微分方程组：
u′ = f ( t , u )
讨论绝对稳定性是困难的，通常考虑它在解u临近的线性化方程：
∂f ( t , u ) ′ (u u ) = (u u ) ∂u u′ = Au
为此我们讨论线性常微分方组： （1.5.1） 其中A为常数矩阵。 设A可对角化：
λ j 是实或复特征值，作变换u=Tv, 则
（1.5.2）
由于 f (t , u ) = μ u ,故方法(1.3.2)中右端是关于un+j(j=0,…,k)的
线性函数。 此时线性k步法(1.3.2)化为k阶线性齐次差分方程
∑α
j =0
k
j
un + j = h ∑ β j f n + j
j =0
j
k
= h∑ β j μ un + j（1.5.3）
j =0
k
或 其中 h =
∑ (α
k
μh 。
j =0
− h β j ) un+ j = 0
（1.5.4）
在实际计算中，每一步计算总会产生舍入误差， 从而得到 的是un的近似值 u n , 其满足
∑ (α
j =0
k
j
− h β j ) un + j = η n
（1.5.5）
其中 η n 为局部舍入误差。 假设 η n ≤ M , 令 en = un − un 显然 en 满足
这个例题说明步长h=0.00125时各点误差很小,而h=0.125 时误差增长很快，它表明零稳定方法如 h取得不合适计算也是 不稳定的，甚至于有的方法尽管它是零稳定的但计算时不管步 长取多么小都不能使它计算稳定,如 Milne方法。 因此,多步法 （1.3.2）还应研究，当步长h固定时计算的稳定性。 所谓计算(数值)的稳定性与(零)稳定性差别是h不趋于零， 是固定的值。 另外这样定义依赖于f ,使对方法的稳定性分析 产生困难。 为解决这一困难, 使判断方法的稳定性不依赖于f , 通常都把方法应用于实验方程。
（1.5.12）
命题5.1 （绝对稳定性的必要条件）一个相容、稳定的多步法 若绝对稳定， 则 μ 的实部 Reμ < 0 , 从而绝对稳定域 D A ⊆ Z − 。
故 1 必为 ρ (λ ) 的单根 ， 设k阶代数方程 ρ (λ ) − h σ (λ ) = 0 的根为：
证明： 因多步法相容、稳定，即
有很多判别法,有著名的Routh-Hurwitz准则, Schur准则 和Miller准则。 下面给出简单而常用的判别法
实系数二次方程 λ2 − bλ − c = 0的根在单位圆内的充分必要条件为:
b < 1− c < 2
实系数三次方程 的充分必要条件为:
λ3 + bλ2 + cλ + d = 0的三个根按模小于1
un + 3 = un + 2 +
的稳定性多项式为：
h (23 f n + 2 − 16 f n +1 + 5 f n ) 12
3 2
1 当h=0.0012时， h = − ，故特征方程为 8 73 2 16 5 3 λ − λ − λ+ =0 96 96 96
h ρ (λ ) − h σ (λ ) = λ − λ − ( 23λ 2 − 16λ + 5 ) , h = −100 h 12
(e
h
− λ1 (h )
) (e
h
− λ2 (h )
)L
(
e h − λk ( h )
)
= O ( h p +1 )
当 h → 0 。故除第一个因子外,其余因子均有正下界, 即 因 eh → 1，
λ j ( h ) − 1 ≥ δ0 > 0
因此
j = 2, 3, L , k （δ 0 是与 h 无关的正数）
( λ < 1 ) 。若存在复数
域DA,使多步法对 ∀h ∈ D A 都绝对稳定, 则称DA为绝对稳定域。
β ) 为绝对稳定区域（区间）。 称DA与实轴的交 (α ，
绝对稳定域 设DA是多步法的绝对稳定域，则稳定多项式确定一由单位圆
λ < 1 到复平面 h ∈ Z 的解析变换：
ρ (λ ) h= σ (λ )