武汉大学：《数值分析》课程教学资源（PPT课件讲稿）第八章 刚性方程组及其数值计算续

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

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

考虑如下线性常微分方程组 y=My 3 y∈R y(0)=(2,1,2), 其中 0.1-4990 M=0-50 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()=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-Kutta方法求解上述初 值问题要求计算直至得到符合精度要求的稳态 解为止我们讨论计算过程可能遇到的问题 一稳定性要求 1|≤278,z=1,2,3 2.78 123|=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( ) 应满足

z()=fy(t,y(t)z(t),0≤t≤7, (1)在y()处的线性化方程 记矩阵J(t)=f(t,y(t)的特征值为 旯(t),i=1,2,…,n 若条件 Rex(t)<0,i=1,2,…,n.0≤1≤T.(3) 称解y()为局部稳定,否则是不稳定的

( ) ( , ( )) ( ), 0 , (2) y z t f t y t z t t T   =   (1)在 y t( ) 处的线性化方程 记矩阵 ( ) ( , ( )) y J t f t y t =  的特征值为 ( ), 1, 2,..., . i  t i n = 若条件 Re ( ) 0, 1, 2,..., . 0 . (3) i  t i n t T  =   称解 y t( ) 为局部稳定,否则是不稳定的

定义:设y(t),t∈[O,7]是方程组(1)的一个解 假定相应 Jacobi矩阵J(t)的特征值满足(3),并且 maxEd(l>>min reM (t), VtE[O,T 则称()在y()附近为刚性方程组刚性比 marEn(t) minRe a(o)

定义:设 y t t T ( ), [0, ]  是方程组(1)的一个解. 假定相应Jacobi矩阵 J t( ) 的特征值满足(3),并且 max Re ( ) Re ( ) , [0, ] i i min i i   t t t T    则称(1)在 y t( ) 附近为 刚性方程组. Re ( ) Re ( ) max min i i i i t t   刚性比

刚性方程组的解 快变部分 慢变部分 般地说,隐型方法比显型方法具有更大 的绝对稳定区域,因此使用隐型方法求解 刚性方程组更为合适 隐型 Runge-Kutta方法 Adams内插方法

刚性方程组的解 快变部分 慢变部分 一般地说,隐型方法比显型方法具有更大 的绝对稳定区域,因此使用隐型方法求解 刚性方程组更为合适. 隐型Runge-Kutta方法 Adams内插方法