第八章 刚性 (stiff)方程组及其数值计算武汉大学数学与统计学院一个化学反应系统中提出的刚性问题的例子
A,B,C是三个化学样本,Robertson反应如下：
7
4
0.04
3 1 0
10
()
()
()
A B s lo w
B B C B v e r y fa s t
B C A C fa s t
我们通过建模可以得到如下方程组
'41 1 2 30.0 4 10y y y y
1 (0) 1y?
' 4 7 22 1 2 3 20.04 10 3 10y y y y y2 (0) 0y?
'3y? 72
23 10 y?
3 (0) 0y?
A：
B：
C：
1 2 3,,Ty y y y?令则方程组可写为：
' ()
(0 ) [1,0,0 ] T
y f y
y
44
32
4 7 4
3 2 2
7
2
0,04 10 10
0,04 10 6 10 10
0 6 10 0
yy
y y y
y
对此 Robertson系统而言,Jacobian矩阵为在平衡点（ 1,0.0000365,0）,Jacobian
矩阵为
0.0 4 0 0.3 65
0.0 4 21 90 0.3 65
0 21 90 0
1 0
2 0,4 0 5
3 2 1 8 9,6
特征值为
，
，
考虑如下线性常微分方程组,
3,,
(0 ) ( 2,1,2 ),T
y M y y R
y
其中
0.1 49.9 0
0 50 0
0 70 30000
M
这里矩阵 M的特征值为
1 2 30,1,5 0,3 0 0 0 0,
上述初值问题的精确解是,
0,1 5 0
1
50
2
5 0 3 0 0 0 0
3
( ),
( ),
( ),
tt
t
tt
y t e e
y t e
y t e e
显然当 t 时解的各个分量 ( ),1,2,3
iy t i?
是指数衰减的,并趋于稳态解
1 2 3(,,) ( 0,0,0 ),y y y?
1 2 3( ),( ),( )y t y t y t
趋于稳态解 (0,0,0 ) 的速度是由因子
0.1te? 决定的,
假如试图利用四级 Runge-Kutta方法求解上述初值问题,要求计算直至得到符合精度要求的稳态解为止,我们讨论计算过程可能遇到的问题,
一,稳定性要求,
4
3
2,7 8,1,2,3,
2,7 8
3 0 0 0 0 1 0,
30000
ihi
h
二,为使解充分接近稳态解只需要,
0,1 0.te
0,1 4
40.
tee
t
而实际上 1t? 后
50 30 00 0,ttee
已经不起作用了 !!!
5
4
40
40 40
4 10
10
t
N
h?
往后的计算我们当然希望使用大步长 !但由于稳定性要求,仍要用小步长,从而耗费了巨大的计算量,并且误差积累的影响也随着计算步数的增加越来越严重,
N是计算步数
410N?
上述例子中,系数矩阵的特征值虽然都是负数,
但绝对值相差非常悬殊,
考虑 n维非线性常微分方程组
( ) (,( ) ),0,( 1 )
( ),[0,],,[0,],n n n
y t f t y t t T
y t T R f T R R
设 ()yt 是 (1)定义在 [0,T]上的解,并满足
0( 0 ),yy?
现用 ( ) ( ) ( )y t y t z t 表示 (1)在 ()yt 附近的解,
则 ()zt 应满足
( ) (,( ) ) ( ),0,( 2 )yz t f t y t z t t T
(1)在
()yt
处的线性化方程记矩阵
( ) (,( ) )yJ t f t y t
的特征值为
( ),1,2,.,,,.i t i n
若条件
R e ( ) 0,1,2,.,,,,0,( 3 )i t i n t T
称解 ()yt 为局部稳定,否则是不稳定的,
定义,设 ( ),[ 0,]y t t T? 是方程组 (1)的一个解,
假定相应 Jacobi矩阵
()Jt
的特征值满足 (3),并且
R e ( ) R e ( ),[ 0,]m a x m i nii
ii
t t t T
则称 (1)在 ()yt 附近为 刚性方程组,
R e ( )
R e ( )
m a x
m i n
i
i
i
i
t
t
刚性比刚性方程组的解快变部分 慢变部分一般地说,隐型方法比显型方法具有更大的绝对稳定区域,因此使用隐型方法求解刚性方程组更为合适,
隐型 Runge-Kutta方法
Adams内插方法
A,B,C是三个化学样本,Robertson反应如下：
7
4
0.04
3 1 0
10
()
()
()
A B s lo w
B B C B v e r y fa s t
B C A C fa s t
我们通过建模可以得到如下方程组
'41 1 2 30.0 4 10y y y y
1 (0) 1y?
' 4 7 22 1 2 3 20.04 10 3 10y y y y y2 (0) 0y?
'3y? 72
23 10 y?
3 (0) 0y?
A：
B：
C：
1 2 3,,Ty y y y?令则方程组可写为：
' ()
(0 ) [1,0,0 ] T
y f y
y
44
32
4 7 4
3 2 2
7
2
0,04 10 10
0,04 10 6 10 10
0 6 10 0
yy
y y y
y
对此 Robertson系统而言,Jacobian矩阵为在平衡点（ 1,0.0000365,0）,Jacobian
矩阵为
0.0 4 0 0.3 65
0.0 4 21 90 0.3 65
0 21 90 0
1 0
2 0,4 0 5
3 2 1 8 9,6
特征值为
，
，
考虑如下线性常微分方程组,
3,,
(0 ) ( 2,1,2 ),T
y M y y R
y
其中
0.1 49.9 0
0 50 0
0 70 30000
M
这里矩阵 M的特征值为
1 2 30,1,5 0,3 0 0 0 0,
上述初值问题的精确解是,
0,1 5 0
1
50
2
5 0 3 0 0 0 0
3
( ),
( ),
( ),
tt
t
tt
y t e e
y t e
y t e e
显然当 t 时解的各个分量 ( ),1,2,3
iy t i?
是指数衰减的,并趋于稳态解
1 2 3(,,) ( 0,0,0 ),y y y?
1 2 3( ),( ),( )y t y t y t
趋于稳态解 (0,0,0 ) 的速度是由因子
0.1te? 决定的,
假如试图利用四级 Runge-Kutta方法求解上述初值问题,要求计算直至得到符合精度要求的稳态解为止,我们讨论计算过程可能遇到的问题,
一,稳定性要求,
4
3
2,7 8,1,2,3,
2,7 8
3 0 0 0 0 1 0,
30000
ihi
h
二,为使解充分接近稳态解只需要,
0,1 0.te
0,1 4
40.
tee
t
而实际上 1t? 后
50 30 00 0,ttee
已经不起作用了 !!!
5
4
40
40 40
4 10
10
t
N
h?
往后的计算我们当然希望使用大步长 !但由于稳定性要求,仍要用小步长,从而耗费了巨大的计算量,并且误差积累的影响也随着计算步数的增加越来越严重,
N是计算步数
410N?
上述例子中,系数矩阵的特征值虽然都是负数,
但绝对值相差非常悬殊,
考虑 n维非线性常微分方程组
( ) (,( ) ),0,( 1 )
( ),[0,],,[0,],n n n
y t f t y t t T
y t T R f T R R
设 ()yt 是 (1)定义在 [0,T]上的解,并满足
0( 0 ),yy?
现用 ( ) ( ) ( )y t y t z t 表示 (1)在 ()yt 附近的解,
则 ()zt 应满足
( ) (,( ) ) ( ),0,( 2 )yz t f t y t z t t T
(1)在
()yt
处的线性化方程记矩阵
( ) (,( ) )yJ t f t y t
的特征值为
( ),1,2,.,,,.i t i n
若条件
R e ( ) 0,1,2,.,,,,0,( 3 )i t i n t T
称解 ()yt 为局部稳定,否则是不稳定的,
定义,设 ( ),[ 0,]y t t T? 是方程组 (1)的一个解,
假定相应 Jacobi矩阵
()Jt
的特征值满足 (3),并且
R e ( ) R e ( ),[ 0,]m a x m i nii
ii
t t t T
则称 (1)在 ()yt 附近为 刚性方程组,
R e ( )
R e ( )
m a x
m i n
i
i
i
i
t
t
刚性比刚性方程组的解快变部分 慢变部分一般地说,隐型方法比显型方法具有更大的绝对稳定区域,因此使用隐型方法求解刚性方程组更为合适,
隐型 Runge-Kutta方法
Adams内插方法
