武汉大学信息与计算科学系：《数值分析》第八章 一阶常微分方程初值问题的数值方法（8-2）常微分方程组

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

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

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

考虑如下线性常微分方程组: 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

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

上述初值问题的精确解是: 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方法求解上述初 值问题,要求计算直至得到符合精度要求的稳态 解为止我们讨论计算过程可能遇到的问题 稳定性要求 内≤278,i=1,2,3 2.78 30000 h40

假如试图利用四级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> 5 ≈ 4×10 h10 而实际上t>1后 N>104 NN,23000经不起作用了 往后的计算我们当然希望使用大步长!但由 于稳定性要求,仍要用小步长从而耗费了巨 大的计算量,并且误差积累的影响也随着计 算步数的增加越来越严重

而实际上 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)定义在0,T上的解,并满足 DOO)=yo 现用y()=y(t)+z()表示()在y()附近的解 则z()应满足

上述例子中,系数矩阵的特征值虽然都 是负数,但绝对值相差非常悬殊. 考虑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(t)=f(t,y(t)z(t),0≤t≤T,(2) (1)在y(t)处的线性化方程 记矩阵J(t)=f(t2y(t)的特征值为 旯(t),i=1,2,……,n. 若条件 Re(t)<0,i=1,2,,n.0≤t≤.(3) 称解y(t)为局部稳定否则是不稳定的

( ) ( , ( )) ( ), 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∈[O,7]是方程组(1)的一个解 假定相应 Jacobi矩阵J(t)的特征值满足(3),并且 max Reλ(> min Re4(),vt∈[0,7 则称()在y()附近为刚性方程组不刚性比 marEn,(t) minren(t)

定义:设 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方法 AdamsI内插方法

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