第2章常微分方程数值解法

1、兰州交通大学数理与软件工程学院 第6章 常微分方程数值解法 兰州交通大学数理与软件工程学院 引言 在科学技术和工程实际问题中，常常需 要求解常微分方程。我们学习过几种类型 的常微分方程的解析解求解方法，在更多 情况下，无法求得解析解。而在多数工程 应用问题中，往往不一定要求解析解，只 需知道解在若干点上的函数值，即求数值 兰州交通大学数理与软件工程学院 ? 常微分方程初值问题 求其数值解，就是计算出解函数()在离 散点 处的近似值 。 ? ? ? ? ? ? ? ? ? ? yf x yaxb y ay , 0 xxxn 12 ? y yyn 12 ,? 兰州交通大学数理与软件工程学院 6.1

2、 初值问题的Euler方法 ,.)2 , 1 , 0()( ,)( 0,.),2 , 1 , 0( )( ),( 00 ? ? ? ? ? ? ? ? ? nyxy yxxy hnnhax yxy yxf dx dy nn nn n 即上的近似值在点去获得解 值解是指通过某种方法假定为常数。该式的数 为步长，一般总记 问题设一阶常微分方程初值 兰州交通大学数理与软件工程学院 初值问题的Euler方法 。算出 出发，逐次公式，它可以从这就是显式的 ）（ 的近似值，则有表示以 于是该式可离散为：代替散化，并用 方法首先将微分算子离为实现这一目标， ., Euler )1 (,.2, 1 ,0),(

3、 )( )(,( )()( , Euler 321 0 1 0 yyy y nyxhfyy xyy xyxf h xyhxy xx nnnn nn nn nn n ? ? ? ? 兰州交通大学数理与软件工程学院 初值问题的Euler方法 。 才能得到步要解函数方程的不同在于，它每算一 方法，它与显式公式或向后这就是隐式的 ）（ 的近似值，则有表示以 于是该式可离散为：代替如果用 1 111 11 01 )2( EulerEuler )2(,.2 , 1 , 0),( )( )(,( )()( , ? ? ? ? ? ? ? n nnnn nn nn nn n y nyxhfyy xyy xyx

4、f h xyhxy xx 兰州交通大学数理与软件工程学院 初值问题的Euler方法 )1( 11 )()1( 1 )( 11 ) 1( 1 )0( 1 111 ,| 3 ,.2 , 1 , 0),(),( 2 ),( ,.2 , 1 , 0),(),( 2 ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? k nn k n k n k nnnnn k n nnnn n nnnnnn yyyy kyxfyxf h yy yxhfyy y nyxfyxf h yy 取时当 ）（ ）（ 时常用以下迭代式：计算 称为梯形公式。 ）（ 平均值的结果，则得如果取以上两式的算术 ? 兰

5、州交通大学数理与软件工程学院 初值问题的Euler方法 收敛。 ）产生的序列时，则由（即 满足常数。如果步长为条件， 满足对变量设函数定理 .)2 , 1 , 0(3 2 , 1 2 0Lipschitz Lipschitzy),(1 . 1 . 6 )( 1 ? ? ? ky L h hL hL yxf k n 兰州交通大学数理与软件工程学院 初值问题的Euler方法 。故有由假设知 ）有）和（证明：由式（ 1 )1( 1 1 )0( 11 1 )( 11 )( 1111 1)(k 1n1n lim,0) 2 (lim: |) 2 ( . | 2 | ),(),(| 2 |yy| 32 ?

6、? ? ? ? ? ? ? ? ? ? ? ? ? ? ? n k n k k k nn k k nn k nnnn yy hL yy hL yy hL yxfyxf h 兰州交通大学数理与软件工程学院 初值问题的Euler方法 。并取 即校正算法构成一类预估迭代一次 一般只由于迭代工作量较大计算对于 )c( 11 )( 11 )( )( 1 1 ),(),( 2 ),( , ,)2( ? ? ? ? ? ? ? ? ? ? ? ? ? nn p nnnnn c n nnn p n n yy yxfyxf h yy yxhfyy y 兰州交通大学数理与软件工程学院 初值问题的Euler方法 )

7、,(,()( 2 ,Euler ,.)2 , 1 , 0(),( ),( )( 2 1 11 12 1 211 nnnnnnnn nn nn nn yxhfyxfyxf h yy nkyhxhfk yxhfk kkfyy ? ? ? ? ? ? ? ? ? ? ? ? ? 亦可写成方法该式称为改进 上式还常写成 兰州交通大学数理与软件工程学院 初值问题的Euler方法 进行比较。精确解 并与法求解法和改进试分别用 设初值问题例 xy y y x y dx dy 21 ,EulerEuler 1)0( 2 1 . 1 . 6 ? ? ? ? ? ? ? ? 兰州交通大学数理与软件工程学院 初值问

8、题的Euler方法 计算结果如下表所示： 法：改进的 法： 上结果，此时计算解：取 ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ,.)2 , 1 , 0() ) 1 . 0(2 ( 1 . 0 ) 2 ( 1 . 0 )( 2 1 Euler ,.)2 , 1 , 0() 2 ( 1 . 0 1 , 0, 1 . 0 1 12 1 211 1 n ky x kyk y x yk kkyy n y x yyyEuler xh n n n n n n nn n n nnn 兰州交通大学数理与软件工程学院 x Euler法y 改进的Euler法y 精确解 0 1.00

9、0000 1.000000 1.000000 0.1 1.000000 1.095909 1.095445 0.2 1.191818 1.184097 1.183216 0.3 1.277438 1.266201 1.264911 0.4 1.358213 1.343360 1.341641 0.5 1.435133 1.416402 1.414214 0.6 1.508966 1.485956 1.483240 0.7 1.580338 1.552514 1.549193 0.8 1.649783 1.616475 1.612452 0.9 1.717779 1.678166 1.67332

10、0 1.0 1.784770 1.737867 1.732051 兰州交通大学数理与软件工程学院 6.1.2 误差概述 有关，称为增量函数。与函数 为而隐式单步法一般形式 显式单步法一般形式为 ),( ),( ),( 111 1 yxf hyxyxhyy hyxhyy nnnnnn nnnn ? ? ? ? ? ? ? 兰州交通大学数理与软件工程学院 误差概述 上的局部截断误差。称为点则 由单步法计算出值计算没有误差，即 的步在点果第上的整体截断误差。如称为在点 则的值或隐式逐步计算，得 出发，由单步法显式从初值定义 11 11 1 1 11111 00 ,)( ,),( )(, )(6.1.

11、1 ? ? ? ? ? ? ? ? nnnn nnn nn nnnnn xyxyT yxyy xnx yxyeyx yxy 兰州交通大学数理与软件工程学院 误差概述 称为为局部截断误差主)(,(则非零项 )()( 阶精度。如果则称该算法具有 误差为 截断如果给果给定的算法的21. 6定义 1 21 1 1 1 ? ? ? ? ? ? ? p nn pp nnn p n hxyxg hOhx,yxgT p )O(hT 兰州交通大学数理与软件工程学院 误差概述 )()( 2 )O()(y 2 Euler式对于 )O(h)( 2 Euler对于显式 则精度相对的越高。 越大，误差阶，一个算法，局部截

12、断一般 4 3 1 3 2 1 3 2 1 hOxy h T hx h T xy h T nn nn nn ? ? ? ? ? ? ? ? ? ? ：对于梯形公式 ：法隐 ：法 说来 兰州交通大学数理与软件工程学院 6.1.3 数值稳定性分析 稳 值不大， 则，则称该算法的舍入误的舍入 误差 。如果算则称该算法是数值稳定 | ，得 值方法又算，如果结果为 计算， 则时有一舍入 误差假设 1 111 数 ，也成绝对稳定 且 通过某种数 实际计算 nn nnn nnn nn yy yy y ? ? ? ? ? 兰州交通大学数理与软件工程学院 数值稳定性分析 ? 定义6.1.3 若某数值算法的绝对稳

13、定性区域 包含h平面上的左半平面Re(h)0，则称 该方法是A ? 隐式Euler法是A稳定的。 兰州交通大学数理与软件工程学院 6.2 Runge-Kutta方法 。这里仍假定 ，使局部截断误差，适当选择参数 可设为方法启发，更一般算式受改进的 )(),()( ,.)2 , 1 , 0( ),( ),( Euler 3 111 21 12 1 22111 nnnnn nn nn nn xyyhOyxyT n kyhxhfk yxhfk kkyy ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? 兰州交通大学数理与软件工程学院 Runge-Kutta方法 )(),(),(),()(

14、)(),(),(),( :Taylor 32 3 1 2 2 hOyxfyxfyxfhxyh hOyxfhkyxfhyxhfk nnynnnnxn nnynnxnn ? ? ? ? 展开式由二元函数 )(),(),( ),()()()( 32 2 2211 hOhyxfyxf yxfxyhxyy nnynn nnxnnn ? ? ? ? ? 于是 兰州交通大学数理与软件工程学院 Runge-Kutta方法 一。自由参数，即解答不唯 程，因此有一个由于四个参数，三个方 展式相比较得与 ? ? ? ? ? ? ? ? ? ? ? ? 2 1 2 1 1 :T aylor 2 2 21 ? ? ?

15、? 兰州交通大学数理与软件工程学院 方法。这是改进的 此时算式为可得取 Euler ),( ),( )( 2 1 , 1, 2 1 , 2 1 ) 1 ( 12 1 211 21 ? ? ? ? ? ? ? ? ? ? ? ? kyhxhfk yxhfk kkyy nn nn nn ? 兰州交通大学数理与软件工程学院 Runge-Kutta方法 .K-R ) 2 1 , 2 1 ( ),( , 2 1 , 1, 0)2( 12 1 21 21 方法这是二阶 此时算式为可得取 ? ? ? ? ? ? ? ? ? ? ? ? kyhxhfk yxhfk kyy nn nn nn ? 兰州交通大学数

16、理与软件工程学院 方法。这也是二阶 （ 又有算式可得取 K-R ) 3 2 , 3 2 ( ),( )3 4 1 , 3 2 , 4 3 , 4 1 )3( 12 1 211 21 ? ? ? ? ? ? ? ? ? ? ? ? kyhxhfk yxhfk kkyy nn nn nn ? 兰州交通大学数理与软件工程学院 6.2.2 四阶Runge-Kutta方法 434241432313212 34324214143 23213133 12122 1 443322111 ,.)2 , 1 , 0( ),( ),( ),( ),( KR ? ? ? ? ? ? ? ? ? ? ? ? ? ? ?

17、 ? ? ? ? ? ? ? ? ，约定 分方程初值问题可设为方法类似，解一阶常微与二阶 n kkkyhxhfk kkyhxhfk kyhxhfk yxhfk kkkkyy nn nn nn nn nn 兰州交通大学数理与软件工程学院 四阶Runge-Kutta方法 稳定。但不是该公式结构整齐 算式求得经典四阶类似以上推导过程，可 A, ),( ) 2 1 , 2 1 ( ) 2 1 , 2 1 ( ),( )22( 6 1 KR 33 23 12 1 43211 ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? kyhxhfk kyhxhfk kyhxhfk yxhfk

18、 kkkkyy nn nn nn nn nn 兰州交通大学数理与软件工程学院 6.2.3 R-K法的稳定性 ? ? ? ? ? ? ? ? ? ? ? ? ? ? n n n n yhhhhk yhhhk yhhk yhk y x y )( 4 1 )( 2 1 )( )( 4 1 )( 2 1 ( )( 2 1 ( d d ,K-R 432 4 32 3 2 2 1 ? ? ? ? ?有应用于实验方程法为例以经典 兰州交通大学数理与软件工程学院 R-K法的稳定性 nn n nn hhh h y hhh h kkkkyy ? ? ? ? ? !4 )( ! 3 )( !2 )( )(1 : !

19、4 )( ! 3 )( !2 )( )(1 )22( 6 1 432 1 432 43211 ? ? ? ? ? 相应的误差为 于是 兰州交通大学数理与软件工程学院 R-K法的稳定性 . , , 0 78 . 2 :, 1| ! 4 )( ! 3 )( ! 2 )( )(1| :|,| 432 1 稳定性必须很小才能保证算法 即限制很大步长的绝对值较大时当因此 可得稳定性区间为负实数时当 于是得稳定性区域为由稳定性要求 h h y f h hhh h nn ? ? ? ? ? ? ? ? ? ? ? ? 兰州交通大学数理与软件工程学院 6.2.5 隐式R-K法 ., ,., A, Euler

20、. .,|, ,K-R 1 1 又能保证数值稳定精度提高 既可以构造高阶隐式方法因此限制较小即对步长定的 稳它们均是但从稳定性分析知迭代求解是它们的缺点 需要多次法和梯形方法计算隐式格式的向后 否则数值不稳定 有较大的限制对步长较大时当但从稳定性分析知 使用方便可以直接计算出法的优点是从显式 h y h y f yy n nn ? ? ? ? 兰州交通大学数理与软件工程学院 隐式R-K法 ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? 22212 12111 22212122 21211111 22111 ),( ),( KR ? ? ? ? ? 约定 方法形式为：设二阶隐式 kk

21、yhxhfk kkyhxhfk kkyy nn nn nn 兰州交通大学数理与软件工程学院 隐式R-K法 ? ? ? ? ? ? ? ? ? ? ? ? ? ? 6 1 )()( 3 1 2 1 1 T aylor),( Taylor),( 22221121221111 2 22 2 11 2211 21 1 21 ? ? ? ? 的同次幂系数得：比较 展式的在再反复带入上式，并与 ，展开求得的做二元函数在点，将 h yxy kyxkk nnn jnn 兰州交通大学数理与软件工程学院 隐式R-K法 ? ? ? ? ? ? ? ? ? ? ? ? ) 3 2 , 3 2 ( ),( )3( 4

22、1 K-R 3 2 0) 1 ( 22 211 211 21 kyhxhfk kkyxhfk kkyy nn nn nn 方法，得二阶三级隐式，令? 兰州交通大学数理与软件工程学院 隐式R-K法 ? ? ? ? ? ? ? ? ? ? ? ? ? ? ) 4 1 ) 6 3 4 1 (,) 6 3 2 1 ( ) 6 3 4 1 ( 4 ,) 6 3 2 1 ( )( 2 1 K-R 6 3 2 1 6 3 2 1 )2( 212 2 1 1 211 21 kkyhxhfk k k yhxhfk kkyy nn nn nn 方法得二阶四级隐式，令? 兰州交通大学数理与软件工程学院 6.3 线形

23、多步法 ? 单步法主要依据yn的信息去计算yn+1。线性多 步法是想依据yn，yn-1，yn-r(r1)的信息去计 算yn+1。 ? 考虑到线性组合较为方便，因此，线性多步法一 般形式可设为 .,0 ).,.,1 , 0 , 1)(,(, ).(. 1 10 0111101 否则称为隐式称为显式时当 其中 ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? rkyxff fhy fffhyyyy knknkn r k knh r k knk rnrnnrnrnnn 兰州交通大学数理与软件工程学院 6.3.1 基于数值积分的方法 ? ? ? ? ? ? ? ? ? ? ? ? 1 0 0

24、 )(,()()( )(,()()( )( ),( 01 1 0 00 x x x x dttytfxyxy xx dttytfxyxy yxy bxayxf dx dy 得取 的等价积分方程形式为 初值问题 兰州交通大学数理与软件工程学院 基于数值积分的方法 公式。这就是著名的 故有只能是近似值由于每一步得到的 略去误差项有形公式对右端定积分采用左矩 有对于一般区间 Euler ,.)2, 1 , 0(),( ,)( )(,()()( ,) 1 ( )(,()()( , 1 1 1 1 1 1n ? ? ? ? ? ? ? ? ? ? nyxhfyy xy xyxhfxyxy dttytfx

25、yxy xx nnnn n nnnn x x nn nn n 兰州交通大学数理与软件工程学院 基于数值积分的方法 有望提高精度。 则值多项式代替被积函数）如果对右端用高次插（ 。这就是熟知的梯形公式 （ 略去误差项有公式对右端定积分采用梯形 ,3 ,.)2 , 1 , 0(),(),( 2 ,)2( 111 ? ? nyxfyxf h yy nnnnnn 兰州交通大学数理与软件工程学院 基于数值积分的方法 ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? 1 0 1 0 1 )1(2 11 0 11 rr11 ) 1,.,1 , 0( ! ) 1).(

26、1( ) 1( ),(),()( )()()( :Newton ),(),.,(),(1 rrjds j jsss ds j s b xxfhbxT xTfbhxyxy fxfxfxr j j nrn rr r r j n j jnn nnnnnn ?其中 后插值公式代入得到 由个点：取 兰州交通大学数理与软件工程学院 基于数值积分的方法 算式时称为四阶显式 则），并逐次取式中的代替用 Adams3 )( 720 251 )( )9375955( 24 )( 8 3 )()51623( 12 )( 12 5 )()3( 2 )( 2 )( ,3210(),( )5(5 1 3211 )4( 4

27、 1211 )3( 3 111 2 11 j ? ? ? ? ? ? ? ? ? ? ? j yhxT ffff h yy y h xTfff h yy y h xTff h yy y h xThfyy jfxyy nnnnnn nnnnn nnnn nnn nnn ? ? ? ? 兰州交通大学数理与软件工程学院 基于数值积分的方法 .Adams )( 720 19 )5199( 24 :, ),(),.,(),( :1 )5(5 2111 1r1r11 算式称为四阶隐式 则可同上求得次插值多项式做 个点若取 ?yh ffff h yy r fxfxfxr nnnnnn nnnnnn ? ?

28、? ? ? 兰州交通大学数理与软件工程学院 基于数值积分的方法 ? Adams预估校正法 ? 预估 ? 校正 ? 并取 )9375955( 24 321 )( 1? ? nnnnn p n ffff h yy 519),(9 24 21 )( 11 )( 1? ? nnn p nnn c n fffyxf h yy )( 11 c nn yy ? ? 兰州交通大学数理与软件工程学院 6.3.2 基于Taylar展开式的方法 ).(. , 3),.,1 , 0 , 1(),(, ).(. 33011331101 10 0111101 ? ? ? ? ? ? ? ? ? ? ? ? nnnnnnn

29、 knknkn r k knh r k knk rnrnnrnrnnn fffhyyyy rrkyxff fhy fffhyyyy ? ? ? 即设中取其中 式在线形多步法的一般形 兰州交通大学数理与软件工程学院 )3 , 2 , 1 , 0 , 1()()( ! 4 )( )( ! 3 )( )( ! 2 )( )()()( )()()(,( )3 , 2 , 1 , 0()()( ! 5 )( )( ! 4 )( )( ! 3 )( )( ! 2 )( )()()( )()( ,Taylor )()(,()(4 5)5( 4 )4( 3 )3( 2 6)5( 5 )4( 4 )3( 3 2

30、 45 0 ? ? ? ? ? ? ? ? ? ? ? khOxy kh xy kh xy kh xykhxy khxyxyxyxf khOxy kh xy kh xy kh xy kh xykhxy khxyxy hhxx xyxyxfxyp n nnnn nknknkn nn nnnn nkn knknknkn 即项和处展到展式在用 及阶精度，将为使算法达到 兰州交通大学数理与软件工程学院 基于Taylar展开式的方法 ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? 110832448116 1271233278 1642294 132 1 ,taylo

31、r)(, ),()(,(),( 3211321 3211321 3211321 32101321 3210 11 1 ? ? ? ? ? 可得到方程组 的系数展式比较的并与展开式则得到 ，并令将其代入 j nn nnnnnnn hxyy xyxyxffxyyy 兰州交通大学数理与软件工程学院 基于Taylar展开式的方法 自由参数。 个方程，因此有四个个参数由于方程组是有 局部截断误差 59 )()( ! 5 )405 8055-243321 ( )( 6)5( 5 3 211321 111 hOxy h yxyT n nnn ? ? ? ? ? ? ? 兰州交通大学数理与软件工程学院 算式和

32、它的余式估计。阶显式这是 所以有 ， 可解得）取（ Adams4 )()( 720 251 )9375955( 24 24 9 24 37 24 59 24 55 1 ,0,01 6)5(5 1 3211 32100 3211 ? ? ? ? ? ? ? ? ? ? ? ? hOxyhT ffff h yy nn nnnnnn ? ? 兰州交通大学数理与软件工程学院 算式和它的余式估计。阶隐式这是 所以有 ， 可解得）取（ Adams4 )()( 720 19 )5199( 24 24 1 24 5 24 19 24 9 1 ,0,02 6)5(5 1 2111 31010 3213 ? ?

33、? ? ? ? ? ? ? ? ? ? hOxyhT ffff h yy nn nnnnnn ? ? 兰州交通大学数理与软件工程学院 ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? )()( 45 14 )22( 3 4 Milne4 )()( 90 1 )49( 3 算式simpson3 6)5(5 1 2131 6)5(5 1 1111 hOxyhT fff h yy hOxyhT fff h yy nn nnnnn nn nnnnn 算式）（ ）（ 其他一些常见的算式： 兰州交通大学数理与软件工程学院 ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ?

34、 ? ? ? 2),( 8 3 )9( 8 1 )22( 3 4 )6( )()( 40 1 )2( 8 3 )9( 8 1 Hamming5 1 )( 112 )( 1 213 )( 1 6)5(5 1 1121 nn p nnnn C n nnnn P n nn nnnnnn ffyxf h yyy fff h yy hOxyhT fff h yyy 校正 预估 校正系统预估 算式）（ 兰州交通大学数理与软件工程学院 6.4 数值解法得收敛性与稳定性 ? 在学习了初值问题 (6.1) 的各种数值解法以 后，为了能够正确地运用这些方法，有必 要简单了解一下关于常微分方程初值问题 数值解法的收

35、敛性及稳定性。收敛性讨论 的是当步长趋于零时，方法的整体截断 误差是否趋于零的问题；稳定性则是讨论 计算过程中的扰动 (舍入误差) 对计算结果 兰州交通大学数理与软件工程学院 数值解法的收敛性 ? 对于给定的初值问题，如果所采用的数值 解法对任一固定的节点 ，当步长 数值解收敛于精确解 ，则称该数值方法 是收敛得。 前面几节所介绍的显式单步法可以统一写成 其中 称为增量函数， 的具体形式 依赖于方程中的f(x,y)以及离散方式。 xaih i ? h ? 0 ? ? y xi ? yyhx y hin iiii? ? 1 11?,? ? x y h, ,?x y h, , 兰州交通大学数理与软

36、件工程学院 ? ? 定理 设初值问题 (6.1) 的数值解计算公 式为 (6.25) ，且满足 （1）局部截断误差 （2）增量函数 关于变量y满足 Lipschitz条件，则单步法的整体误差是P阶 的，即当 时， 单步法收敛。 ? eO h i p ? ? ? 1 1 ? x y h, , ? ? ? ? E i hy xyO h ii p ,? p ? 1 兰州交通大学数理与软件工程学院 数值解法的稳定性 ? 一个收敛的数值解法，截断误差的影响随 步长的减小而减小。但另一方面，舍入 误差的影响会随步长的减小而增大。在 使用某种数值方法计算的过程中，如果某 步长产生的舍入误差以后不能逐步减弱，

37、 累积起来势必给结果造成难以估量的影响， 这样的数值方法就不宜采用。 兰州交通大学数理与软件工程学院 ? 如果某种数值方法，在节点 处的 值有大小 为 的扰动 (舍入误差) ，而在其后的各节点 值 的扰动都不超过，则称该数值方法是 (绝 对)稳定的。 ? 各种数值方法的稳定性，依赖于算法过程以及方 程的形式，这取决于方程右端的表达式，给讨论 数值方法的稳定性带来困难。为简便起见，我们 以一个简单的微分方程为例，来说明讨论数值解 法稳定性的过程。以下针对微分方程 进行 讨论，其中为常数 (可以是复常数) 。 xiy i ? xj j ? 1 yj y j dy dx y? 兰州交通大学数理与软件

38、工程学院 Euler法显式格式的稳定区域 ? 兰州交通大学数理与软件工程学院 ? 当 时，第i步运算产生的扰动，在 以后的计算中逐步减弱。从而方程 (6.27) 的Euler法显式格式的稳定区域为 。 在复平面上表示以 (-1，0) 为圆心的单位圆内 部。特别当为负实数 (0) 时，稳定域 为 ，计算时步长应满足 11?H 11?H ?20?h 0 2 ?h ? 兰州交通大学数理与软件工程学院 Euler法隐式格式的稳定区域 ? 设在节点 处值 有扰动 ，由此所引起 的 的扰动为 ， 从而 于是 一般地 ? yyhyyHyin iiiii? ? 111 11?,? xiy i? i yi?1

39、? i?1 ? yyH yin iiiiii? ? 1111 011?, ,? ? iii H ? ? 11 ? iii H ? ? 11 ? ? ? ij i j H jni ? ? ? ? 1 12 , ,? 兰州交通大学数理与软件工程学院 ? 从而得方程 (6.27) 的Euler法隐式格式的稳 定区域 ，在复平面上表示以（0， 1）为圆心地单位圆地外部，它较显式格式 地稳定区域大得多。 当为负实数 (0) 时，对任意的 恒成立。这种对 步长选取无限制的稳定区域称为无条件稳 定区域。 11?H hh?0 11,? 兰州交通大学数理与软件工程学院 梯形公式的稳定区域 ? 设在点 处值 有扰

40、动 ，由此引起的 的扰动为 ，类似于前面地讨论可得 由此可得方程的梯形算法的稳定区域为 它表示左半复平面。当为负实数 (0) 时， 为无条件稳定区域。 ? yyh yyin iiii? ? 11 1 2 011?, ,? xiyi ? i yi?1 ? i?1 ? i j i H H jni ? ? ? ? ? ? ? ? ? ? 1 2 2 12 , ,? 22?HH 兰州交通大学数理与软件工程学院 ? 类似地，可以分析RungeKutta法、 Aadms法的稳定区域，由于推导过程复杂， 故从略。在实际计算中，当方程中的表达 式比较复杂时，分析一个解法的稳定性是 困难的，人们常常通过数值试验

41、进行判断。 经验告诉我们：一个不稳定的算法，计算 结果变化急剧，解的误差往往以指数级增 大。若解出现这种不稳定现象，减小步长 再算，如仍不正常，改用其它数值方法。 兰州交通大学数理与软件工程学院 6.5 微分方程组及高阶微分方程的数值解法 算式为四阶方法求解如果用经典的 时 设一阶方程组 ,K-R , ),( d d ),( d d 2021010 212 2 211 1 ? ? ? ? ? ? ? ? ? ? ? ? yyyytt yytf t y yytf t y 兰州交通大学数理与软件工程学院 一阶常微分方程组的R-K方法 ,.)2 , 1 , 0( ),(),( ) 2 , 2 , 2

42、 (), 2 , 2 , 2 ( ) 2 , 2 , 2 (), 2 , 2 , 2 ( ),(),( )22( 6 1 )22( 6 1 232131221232131114 22 2 12 1223 22 2 12 1113 21 2 11 1222 21 2 11 1112 2122121111 24232221212? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? n kykyhthfkkykyhthfk k y k y h thfk k y k y h thfk k y k y h thfk k y k y h thfk

43、 yythfkyythfk kkkkyy kkkkyy nnnnnn nnnnnn nnnnnn nnnnnn nn nn ， 兰州交通大学数理与软件工程学院 。步求出 重复上述过程可逐， 的值；再以，时的，再计算出， ，代入上式，先计算 开始，将计算过程：从 ,.)2 , 1 , 0(, 1 ,0 1211 2121111 211112414 2313122111202 1010 ? ? ? ? ? ? nyy yyyyttn yyttkk kkkkkyy yyttn nn 兰州交通大学数理与软件工程学院 一阶常微分方程组的R-K方法 ? ? ? ? ? ? ? ? ? ? ? ? ? ?

44、? ),.,2 , 1()( ),.,( d d . ),.,( d d ),.,( d d 00 21 212 2 211 1 miyty yyytf t y yyytf t y yyytf t y m ii mm m m m 值问题设一阶常微分方程组初 个方程的情形。到将两个方程的情形推广 兰州交通大学数理与软件工程学院 一阶常微分方程组的R-K方法 :K-R ,),.,( )( ),( d d : ,),.,(,),.,( T 020100 00 T 21 T 21 的向量形式是方法 于是经典的四阶，这里 向量形式为 则方程的若记 btayyy t t t fffyyy m mm y y

45、y yf y fy ? ? ? ? ? ? ? ? ? 兰州交通大学数理与软件工程学院 一阶常微分方程组的R-K方法 ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ,.)2 , 1 , 0 ),( ) 2 , 2 ( ) 2 , 2 ( ),( )22( 6 1 34 2 3 1 2 1 43211 n hth h th h th th kyfk k yfk k yfk yfk kkkkyy nn nn nn nn nn （ , ),( ),( ),( ),( . 1 12 11 2 1 1 2 1 ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? mnnm nn nn nn mm n n n k k k thf thf thf th y y y y y y yfky ? ， 这里向量 兰州交通大学数理与软件工程学院 , ) 2 , 2 ( ) 2 , 2 ( ) 2 , 2 ( ) 2 , 2 ( , ) 2 , 2 ( ) 2 , 2 ( ) 2 , 2 ( ) 2 , 2 ( 3 32 31 2 2 2 2 1 2 3 2 22 21 1 1 2 1