武汉大学数学与统计学院：《数值分析》第八章 常微分方程的数值方法（8.3）stiff systems

第八章刚性(sti方程组及其数值 计算 武汉大学数学与统计学院

第八章 刚性(stiff)方程组及其数值 计算 武汉大学数学与统计学院

个化学反应系统中提出的刚性 可题的例子 A,B,C是三个化学样本, Robertson反应如下 A- >b(slow) B+B-10>C+B(very fast B+C->A+C(fast)

一个化学反应系统中提出的刚性 问题的例子 A，B，C是三个化学样本，Robertson反应如下： 7 4 0.04 3 10 10 ( ) ( ) ( ) A B slow B B C B very fast B C A C fast  ⎯⎯→ + ⎯⎯⎯→ + + ⎯⎯→ +

我们通过建模可以得到如下方程组 :y1=-0.04y1+10y2y3 y(0) B :y2=0.04y1-104y2y3-3·10y2y2(0)=0 C: y3 3.107y2y(0)=0

我们通过建模可以得到如下方程组 ' 4 1 1 2 3 y y y y = − + 0.04 10 1 y (0) 1 = ' 4 7 2 2 1 2 3 2 y y y y y = − −  0.04 10 3 10 2 y (0) 0 = ' 3 y = 7 2 2 3 10  y 3 y (0) 0 = A： B： C：

y=(y,y2,y3 则方程组可写为: y=f() y(0)=[,0,0] 对此 Robertson系统而言, Jacobian矩阵为 0.04 104y 10 0.04-104y3-6:107y2-10y2 6107y2

( 1 2 3 , , ) T 令y y y y = 则方程组可写为： ' ( ) (0) [1,0,0]T y f y y = = 4 4 3 2 4 7 4 3 2 2 7 2 0.04 10 10 0.04 10 6 10 10 0 6 10 0 y y y y y y   −   − −  −        对此Robertson系统而言，Jacobian矩阵为

在平衡点(1,0.000065,0), Jacobian 矩阵为 0.040 0.365 0.04-2190-0.365 02190 0 特征值为21=0 =-0.405 13=-21896

在平衡点（1，0.0000365，0），Jacobian 矩阵为 0.04 0 0.365 0.04 2190 0.365 0 2190 0   −   − −       1  = 0 2  = −0.405 3  = −2189.6 特征值为 ，

考虑如下线性常微分方程组: y=My 3 y∈R y(0)=(2,1,2), 其中 0.1-4990 M=0-500 070-30000 这里矩阵M的特征值为 1=-0.1,2=-50,3=-30000

考虑如下线性常微分方程组: 3 , , (0) (2,1, 2) , T y My y R y  =  = 其中 0.1 49.9 0 0 50 0 0 70 30000 M   − −   = −     −   这里矩阵M的特征值为 1 2 3    = − = − = − 0.1, 50, 30000

上述初值问颎的精确解是: y1(t)=e01+c-501 2(4)=e-507 y3(t)=e-30+ 30000t 显然当t→+∞时解的各个分量y(1),i=1,2,3 是指数衰减的并趋于稳态解(1,y2,y3)=(0,0,0) 1(D),y2(D),y3(t)趋于稳态解(0,0,0)的速度是 由因子e1决定的

上述初值问题的精确解是: 0.1 50 1 50 2 50 30000 3 ( ) , ( ) , ( ) . t t t t t y t e e y t e y t e e − − − − − = + = = + 显然当 t → + 时解的各个分量 ( ), 1, 2,3 i y t i = 是指数衰减的,并趋于稳态解 1 2 3 ( , , ) (0,0,0). y y y = 1 2 3 y t y t y t ( ), ( ), ( ) 趋于稳态解 (0, 0, 0) 的速度是 由因子 0.1t e − 决定的

假如试图利用四级 Runge-Kutt方法求解上述 初值问题要求计算直至得到符合精度要求的稳 态解为止我们讨论计算过程可能遇到的问题 稳定性要求: 1|≤278,z=1,2,3 2.78 |=30000∴h< 10 30000 为使解充分接近稳态解只需要: 0.1t )0 →t4O

假如试图利用四级Runge-Kutta方法求解上述 初值问题,要求计算直至得到符合精度要求的稳 态解为止.我们讨论计算过程可能遇到的问题: 一.稳定性要求: 4 3 2.78, 1, 2,3. 2.78 30000 10 . 30000 h i i h   −  = =    二.为使解充分接近稳态解只需要: 0.1 0. t e − → 0.1 4 40. t e e t − −   

t>40 N是计算步数 4040 →N> 4×10 h10 而实际上t>1后 N>1O4 e30,e300u经不起作用了!! 往后的计算我们当然希望使用大步长!但由 于稳定性要求仍要用小步长.从而耗费了巨 大的计算量并且误差积累的影响也随着计 算步数的增加越来越严重

而实际上 t 1 后 50 30000 , t t e e − − 已经不起作用了!!! 5 4 40 40 40 4 10 10 t N h −     =  往后的计算我们当然希望使用大步长!但由 于稳定性要求,仍要用小步长.从而耗费了巨 大的计算量,并且误差积累的影响也随着计 算步数的增加越来越严重. N是计算步数 4 N 10

上述例子中系数矩阵的特征值虽然都是负数, 但绝对值相差非常悬殊. 考虑n维非线性常微分方程组 y(t)=f(t,y(t),0≤t≤T, y(t):[0,7]→>R,f:[0,7×R”→R 设y(t)是(1)定义在,T上的解,并满足 V(O)=o 现y(t)=y()+z(t)表示(1)在y(t)附近的解 则z(t)应满足

上述例子中,系数矩阵的特征值虽然都是负数, 但绝对值相差非常悬殊. 考虑n维非线性常微分方程组 ( ) ( , ( )), 0 , (1) ( ) :[0, ] , :[0, ] . n n n y t f t y t t T y t T R f T R R  =   →  → 设 y t( ) 是(1)定义在[0,T]上的解,并满足 0 y y (0) . = 现用 y t y t z t ( ) ( ) ( ) = + 表示(1)在 y t( ) 附近的解, 则 z t( ) 应满足