Runge-Kutta法

常微分方程初值问题数值解法 第1章 §5 单步法 Runge-Kutta法 §5 单步法 Runge-Kutta法 Euler是最简单的单步法。单步法不需要附加初值，所需的 存储量小，改变步长灵活，但线性单步法的阶最高为2。本节将 介绍非线性（关于f ）高阶单步法，重点是Runge-Kutta法。 5.1 Taylor展开法 设初值问题 ),( utfu =′ 00 )( utu = 的解充分光滑，将 u(t) 在t0处用Taylor公式展开： ( )10)(0200 )(!)(!2)()()( ++++′′+′+= pp p hOtu p htuhtuhtutu L 其中 ,)( 00 utu = (5.1) ),,())(,()( 00000 utftutftu ==′ 0 0( ) t t dfu t dt = ′′ = = 0 )( 0 )3( ttdt df dt dtu = ⎥⎦ ⎤⎢⎣ ⎡= ),(),(),(),(2),( 00002000000 utfutfutfutfutf uututt ++= ( )200000000 ),(),(),(),( utfutfutfutf uut ++ LL (5.2) 0 0( , )tf t u 0 0 0 0( , ) ( , )uf t u f t u+ 令 舍去余项 1 1 1 1 ( , ( ); ) ( , ( )) ! j jp j j h dt u t h f t u t j dt ϕ − − − = = ∑ 则可将（7-5）改写成为 (5.3) )());(,()()( 10000 ++=−+ phOhtuthtuhtu ϕ ，则得)( 1+phO 1 0 0 0( , ; )u u h t u hφ− = 。 一般而言，若已知un，则 1 ( , ; )n n n nu u h t u hϕ+ = + 这是一个单步法，局部截断误差为O(h p+1)， （5.3），可知f关于 f 非线性。 由（5.2） 当 p=1时，它是Euler法。 由于计算f(tn,un;h)的工作量太大，一般不直接用Taylor 展开法做数值计算，但可用它计算附加值。 0,1, 2, ,n = L (5.4) (5.5) 用4阶Taylor公式求如下初值问题： 分别取h=0.25,0.125并比较结果。 ( ) , 2 t uu t −′ = (0) 1u = 30 ≤< t 首先， ( ) , 2 t uu t −′ = ( ) 2 , 2 4 d t u t uu t dt − − +⎛ ⎞′′ = =⎜ ⎟⎝ ⎠ ( ) 2 2 , 4 8 d t u t uu t dt − + − + −⎛ ⎞′′′ = =⎜ ⎟⎝ ⎠ ( ) ( )4 2 2 , 8 16 d t u t uu t dt − + − + −⎛ ⎞= =⎜ ⎟⎝ ⎠ 要求出u1,必须在点(t0,u0)=(0,1)处的上述各阶导数值，经计算得： 例 ( ) 0.0 1.00 0.5, 2 u −′ = = − ( ) 2 0.0 1.00 0.75,4u − +′′ = = ( ) 2 0.0 1.00 0.375, 8 u − + −′′′ = = − ( ) ( )4 2.0 0.0 1.00 0.1875, 16 u + −= = 然后将导数和h=0.25代入（5.5）式，计算出： ( )2 31 0 0.5 0.75 0.375 0.18752 6 24 h h hu u h ⎛ ⎞= + − + × + × − + ×⎜ ⎟⎝ ⎠ 0.75 0.375 0.18751.0 0.25 0.5 0.25 0.25 0.25 2 6 24 ⎛ ⎞⎛ − ⎞⎛ ⎞= + × − + × + × + ×⎜ ⎟⎜ ⎟⎜ ⎟⎝ ⎠⎝ ⎠⎝ ⎠ 8974915.0= 计算出解点为(t1, u1)=(0.25, 0.8974915)。 要计算u2必须在点(t1,u1)=(0.25,0.8974915)处的上述各阶导数值， 计算量非常大，手算十分枯燥，经计算得： ( ) 0.25 0.89749150.25 0.3237458, 2 u −′ = = − ( ) 2 0.25 0.89749150.25 0.6618729, 4 u − +′′ = = ( ) 2 0.25 0.89749150.25 0.3309364, 8 u − + −′′′ = = − ( ) ( )4 2.0 0.25 0.89749150.25 0.1654682, 16 u + −= = 然后将导数和h=0.25代入（5.5）式，计算出： ( )2 32 1 0.3237458 0.6618729 0.3309364 0.16546822 6 24 h h hu u h ⎛ ⎞= + + × + × + ×⎜ ⎟⎝ ⎠ - - 0.6618729 0.3309364 0.16546820.8974915 0.25 0.3237458 0.25 0.25 0.25 2 6 24 ⎛ ⎞⎛ − ⎞⎛ ⎞= + × − + × + × + ×⎜ ⎟⎜ ⎟⎜ ⎟⎝ ⎠⎝ ⎠⎝ ⎠ 0.8364037= 计算出解点为(t2, u2)=(0.50, 0.8974915)。 经大量的计算后，得到用4阶Taylor公式求如上初值问题的 如下的数值结果，见下 表 关于同志近三年现实表现材料材料类招标技术评分表图表与交易pdf视力表打印pdf用图表说话 pdf 。 h=0.25 h=0.125 精确解 un un u(tn) 0.000 1.000 1.000 1.000 0.250 0.8974915 0.8974908 0.8974917 0.500 0.8364037 0.8364024 0.8364023 0.750 0.8118696 0.8118679 0.8118678 1.000 0.8195940 0.8195921 0.8195920 1.500 0.9171021 0.9170998 0.9170997 2.000 1.1036408 1.1036385 1.1036383 2.500 1.3595168 1.3595145 1.3595144 3.000 1.6693928 1.6693906 1.6693905 tn 5.2 单步法的收敛性 其中p≥1是使(5.7)成立的最大整数，则称算法 Taylor展开法是p阶单步法F( t,u,h)由(5.3)定义，Euler法 是一阶单步法，相应的F=F(t, u(t))。 将初值问题写成积分形式 如果有某已确定的函数F(t, u, h)（通过某种离散化），使初值 问题的任意解u(t) 满足 ( ) ( ) ( )( ) τττ duftuhtu ht t∫ +=−+ , ( ) 00 utu = (5.6) ( ) ( ) ( )( ) ( )1, , ,pu t h u t h t u t h O hφ ++ − = + (5.7) h Tn ,,1,0 L=( ) ( ) ( ), , ,n nu t h u t h t u hφ+ − = (5.8) 为p阶单步法。 注意，由积分中值定理可得： ( )( ) ( 则可知F( t, u(t),h)于h=0连续。 ) ( ) ( ), , pu t h u tt u t h O hhφ + −= + ( )( ) ( )1 , t h p t f u d O h h τ τ τ+= +∫ 故： ( )( ) ( )( ) 0 lim , , , h t u t h f t u tφ→ = ( )( ) ( )phOuf += ξξ , ( )( ) ( )( ), ,0 ,t u t f t u tφ =义 ，定 反之，若F( t, u(t),0)=f(t, u(t), h)于h=0连续，且单步法（5.8）的 局部截段误差Rh=o(h)(h→0)，则 ( )( ) ( )1 , 1t h t f u d o h τ τ τ+ +∫( )( ), ,t u t hφ = ( )( ) ( )( )0 , ,0 ,h t u t f t uφ→ =令 则知, 定义为 定义5.1 单步法（5.8）相容， t 。所以可将单步法的相容性 ( )( ), ,t u t hφ如果 于h=0连续，且 ( )( ) ( )( ), ,0 ,t u t f t u tφ = 。 ( )1 , ,n n n nv v h t v hφ+ = + ( ) ( )1 , , , ,n n n n n ne e h t u h t v hφ φ+ = + −⎡ ⎤⎣ ⎦ ( )1 1n n n ne e hL e hL e+ ≤ + = + ( )( )01n h T t+ ≤ − 与（5.8）式相减，并记en=un-vn,得： 事实上，设由初值v0算出的（5.8）的近似值{vn},则 关于u满足Lipschitz条件（Lipschitz常数为L）,单步法（5.8）稳定 则 所以（5.8）是稳定的。 定理5.1 设F(t,u,h)(t0≤t≤T，0≤h≤h0，u∈（-∞，∞）） ( ) ( )01 0 0e e L T tL n h e e−+≤ ≤ ≤L 定理5.2 相容，则当h→0时，它的数值解un收敛到初值问题的精确解u(tn) ， 只要t0+nh→t, 初值u0→u(t0)。若更有（5.8）是p阶单步法，则还有 误差估计式（5.9）。 由于F关于u满足Lipschtz条件，则 1 1 p n n nLh chε ε ε ++ ≤ + + 两边同时关于n求和，得 设F(t,u,h)满足定理5.1的条件，又单步法（5.8） 证明 不妨设（5.8）是p阶单步法，则初值问题的解u(t)满足 ( ) ( ) ( )( ) ( )11 , , ,pn n n nu t u t h t u t h O hφ ++ = + + 与（5.8）式相减，并记en=u(tn)-un,得： ( )( ) ( ) ( )11 , , , , pn n n n n nh t u t h t u h O hε ε φ φ ++ ⎡ ⎤= + − +⎣ ⎦ 1 1 p n n nLh chε ε ε ++ − ≤ + 再利用Gronwall不等式（1.25），就得到估计式： 特别若取u0=u(t0)，则e0=0，整体截断误差的阶为O(hp)。 nε ≤ ( ) 10 0 n p k k cTh Lhε ε− = + + ∑ ( ) 01TL pn e Lh cThε ε⎡ ⎤≤ + +⎣ ⎦ （5.9） 通过观察我们发现显式Euler法和隐Euler法各用到了u(t) 在[t,t+h]上的一个一阶导数值，它们都是一阶 方法 快递客服问题件处理详细方法山木方法pdf计算方法pdf华与华方法下载八字理论方法下载 。梯形法和 改进的Euler法用到了u(t)在[t,t+h]上的两个一阶导数值，它 们都是二阶方法。而Runge-Kutta型方法是用u(t)在[t,t+h]上 的 f 在一些点的值非线性表示f(t,u(t),h) ，使单步法的局部 截断误差的阶和Taylor展开法相等。 5.3 Runge-Kutta法 Taylor展开法，用在同一点(tn,un)的高阶导数表示f(t, u(t), h)， 这不便于计算。 ii m i nn kchuu ∑ = + += 1 1 ),( 1 ∑ = ++= m j jijnii kbhuhatfk ,,,1 mi L= 0 ,nt t nh= + ( )0,1, , ,n = L Runge-Kutta方法是求解初值问题（3.1）的一类重要的 经典算法，显式Runge-Kutta方法的绝对稳定区域是有限区。 不适用于求解刚性方程，1964年Butcher首先提出了隐式的 Runge-Kutta法，可用于求解刚性方程。 其中 h为积分步长，un为u(tn)的 近似值， B=(bij)(i,j=1,…,m)称为Runge-Kutta矩阵，其中在 权。 Runge-Kutta方法的一般结构 上面公式中未出现的元素定义为零，而a=(a1,a2,…,am)T和 c=(c1,c2,…,cm)T分别称为Runge-Kutta节点和Runge-Kutta 且这些系数要满足以下条件： 1 0, 1 , m i i i c c = ≥ =∑ miabm j iij ,,2,1, 1 L==∑ = 通常可将这些系数排列成以下形式，成为Butcher数表示法： 在公式中要用到m个f(t,u)的值，故称为m级Runge-Kutta法。 法，此时k1=f(tn,un)，则由 计算公式 六西格玛计算公式下载结构力学静力计算公式下载重复性计算公式下载六西格玛计算公式下载年假计算公式 可逐个递推出 k2,…,km 如果j≥i时， bij=0，即为严格下三角矩阵，就是显式R-K方 如果j＞i时， bij=0，而对角元素bii不全为零，这时公式称为 半隐式Runge-Kutta方法，这时ki可表示为： Ba Tc ma a M 1 mmm m bb bb L MOM L 1 111 mcc L1 这是关于ki的非线性方程式，因此每一步要解m个非线性方程式， 若bii均相等，则成为对角隐式Runge-Kutta方法简称DIRK方法， 此时由上式求ki的系数矩阵相同，故每步只求一个非线性方程 一般隐式Runge-Kutta方法的右端包含全部k1,k2,…,km,故 为得到Runge-Kutta公式需要确定公式的系数a,B,c。下面 式的解。 求k1,k2,…,km需要解一个m×m阶的非线性方程组。 我们用Taylor展开思想来构造高阶显式Runge-Kutta公式。 ),( 1 1 iii m j jijnii khbkbhuhatfk +++= ∑− = ,,,1 mi L= 先引进若干记号，首先[t,t+h]取上的m个点: ,321 htttttt m +≤≤≤≤≤= L 令 ,,,2,1 mihathatt iii L=+=+= 其中a,B与h无关, 1 1 , 2, , , i ij i j b a i m − = = =∑ L ∑ = =≥ m i ii cc 1 1,0 ( )20, , , ,Tma a= La ⎟⎟ ⎟⎟ ⎟ ⎠ ⎞ ⎜⎜ ⎜⎜ ⎜ ⎝ ⎛ = −121 3231 21 mmmm bbb bb b L LLB 此时 B为严格三角矩阵。 （5.12） （5.13） 假设三组系数和已给定，则显式Runge-Kutta法的计算过程 L,1,0),,,(1 =+=+ nhuthuu nnnn φ （5.14） 其中 ∑ = = m i ii kchtut 1 )),(,(φ 如下： 1 ( , ),k f t u= 2 2 21 1( , ( ) ) ,k f t ha u t hb k= + + 3 3 31 1 32 2( , ( ) ( )),k f t ha u t h b k b k= + + + 31 32 3,b b a+ = 1 1 ( , ( ) ) , m m m mj j j k f t ha u t h b k − = = + + ∑ L L 21 2 ,b a= 1 1 m mj m j b a − = =∑ （5.15） （5.16） ( ) ( )111 +++ =− pnn hOutu 假定u(tn)=un ，由显式Runge-Kutta公式得到的数值解un+1 与方程的精确解u(tn+1)之差，如果有 系数{ai},{bij}和{ci}将按如下原则确定： 到（5.15）式中，使hl(l=0,1,…,p-1)的系数和Taylor展式 同次幂的系数相等。 现在推导一些常用的计算 方案 气瓶 现场处置方案 .pdf气瓶 现场处置方案 .doc见习基地管理方案.doc关于群访事件的化解方案建筑工地扬尘治理专项方案下载 ，特别地，给出m=3显式 首先将u(t+h)在 t 处展开到 h 的三次幂,即： 3 ( ) 4 1 ( ) ( ) ( ) ( ) ! l l l hu t h u t u t O h l= + = + +∑ (5.17) 其中 (5.18) 如此得到的算法（5.14）称为m级p阶Runge-Kutta法。 Runge-Kutta法的推导。 ( ) ( , ( ), ),u t h t u t hϕ= + % 2 31 1 ˆ( , , ) ( ) ( ) 2 6 u t u h f hf h ff f O hϕ = + + + +% %% t uf f ff= +% 2ˆ 2tt tu uuf f ff f f= + + 将u(tn+1)在tn处Taylor展开，再将ki关于h展开，代入 其次，由二元函数f在点(t,u(t))处的Taylor展开式可得： ,))(,(1 ftutfk == ),( 1222 khauhatfk ++= ( ) )(2 2 1)( 3211 2 2 2 12 hOfkfkfahfkfhaf uututtut ++++++= )(ˆ 2 1~ 32 2 2 2 hOfahfhaf +++= )()ˆ 2 1~(~ 323322 2 33 hOfaffbahfhafk u ++++= 于是，代入（7-21）中，即 ( )3322113 1 ),,( kckckckchut i ii ++== ∑ = φ ( ) 21 2 2 1 31 ˆ2t uc f c f ha f k f a f h⎛ ⎞= + + + +⎜ ⎟⎝ ⎠ 2 2 3 3 3 2 32 3 1 ˆ ( ) 2u c f a f h a b f f a f h O h⎛ ⎞⎛ ⎞+ + + + +⎜ ⎟⎜ ⎟⎝ ⎠⎝ ⎠ % % ( ) ( ) +++++==∑ = fcacahfccckchut i ii ~),,( 3322321 3 1 φ [ ] )(ˆ)(~2 2 1 3 3 2 32 2 23322 2 hOfcacaffcbah u ++++ 由已知 (5.19) (3) 21 1( , , ) ( ) ( ) ( ) 2 6 t u h u t u t h u t hϕ ′ ′′= + +% 其中 ,t uf f ff= +% 2ˆ 2tt tu uuf f ff f f= + + 。 2 31 1 ˆ( ) ( ) 2 6 u f hf h ff f O h= + + + +% % 并合并 hl(l=0,1,2)的 比较 ),,( hutφ 和 ),,(~ hutφ 的同次幂系数，可得以下具体方案： （一）m=1 此时 c1=c2=0, F(t,u,h)=c1 f ,比较h的零次幂，知 ,),,( fhut =φ 方法（5.14）为一级一阶Runge-Kutta法，实际上为Euler法。 （二）m=2,此时 c3=0 ( ) )(ˆ 2 1~),,( 32 2 2 2 2221 hOfcahfchafcchut ++++=φ 与 ),,(~ hutφ 比较1,h的系数，则 2 11 2221 ==+ cacc 它有无穷多组解，从而有无穷多个二级二阶方法。 三个常见的方法是： , 2 1,1,0 221 === acc（1） 称为中点法。 其Butcher表为： 此时 1 2 10 1 2 0 1 2 ,n nu u hk+ = + 1 ( , ),n nk f t u= 2 1 1 1, 2 2n n k f t h u hk⎛ ⎞= + +⎜ ⎟⎝ ⎠ 1 2 2 1 , 1 , 2 c c a= = =（2） 称为改进的Euler法。 其Butcher表为 此时 1 2 0 ( )1 1 2 ,2n n hu u k k+ = + + 1 ( , ),n nk f t u= 2 1( , )n nk f t h u hk= + + 1 1 2 1 1 2 2 1 3 2, , , 4 4 3 c c a= = =（3） 其Butcher表为： 此时 3 4 0 ( )1 1 23 ,4n n hu u k k+ = + + 1 ( , ),n nk f t u= 2 1 2 2, 3 3n n k f t h u hk⎛ ⎞= + +⎜ ⎟⎝ ⎠ 2 3 1 4 1 （三）m=3 比较（5.18）和（5.19），令 1,h,h2 的系数 相等，并注意 ff ˆ,~ 的任意性，得 2 1,1 3322321 =+=++ cacaccc . 6 1, 3 1 33223 2 32 2 2 ==+ cbacaca 四个方程不能完全确定六个系数，因此这是含两个参数的三级 三阶方法类。 常见方案有： 1 2 3 1 3, 0, , 4 4 c c c= = = 2 3 321 2 2, ,3 3 3a a b= = = 。 (1) Heun三阶方法。 此时取 Butcher表为： 3 4 0 ( )1 1 33 ,4n n hu u k k+ = + + 1 ( , ),n nk f t u= 2 1 1 1, 3 3n n k f t h u hk⎛ ⎞= + +⎜ ⎟⎝ ⎠ (5.20) 1 3 1 4 3 2 2 2, 3 3n n k f t h u hk⎛ ⎞= + +⎜ ⎟⎝ ⎠ 2 3 0 0 2 3 1 3 （2）Kutta三阶方法， 1 2 3 1 2 1, , , 6 3 6 c c c= = = 2 3 321 , 1, 22a a b= = = 。 方法为: ),4(6 1 3211 kkkhuu nn +++=+ ),,(1 nn utfk = ), 2 1, 2 1( 12 hkuhtfk nn ++= ).2,( 213 hkhkuhtfk nn +−+= (5.21) 此时 Butcher表为： 4 30 4 1 3 2 3 1 0 3 20 3 1 （3）Nystrom三阶方法， 1 2 3 1 3 3, , , 4 8 8 c c c= = = 2 3 322 2 2, ,3 3 3a a b= = = 。 方法为: 1 1 2 33 3 ,4 2 2n n hu u k k k+ ⎛ ⎞= + + +⎜ ⎟⎝ ⎠ ),,(1 nn utfk = 2 1 2 2, , 3 3n n k f t h u hk⎛ ⎞= + +⎜ ⎟⎝ ⎠ 3 2 2 2, 3 3n n k f t h u hk⎛ ⎞= + +⎜ ⎟⎝ ⎠。 此时 Butcher表为： 4 30 4 1 3 2 3 1 0 3 20 3 1 （四）m=4将（5.18）和（5.19）展开到h3,比较 )3,2,1,0( =ihi 的系数，则含13个待定系数的11个方程，由此得到含两个参数的 四级四阶Runge-Kutta方法类，其中最常用的有以下两个方法： ),22( 6 43211 kkkkhuu nn ++++=+ ),,(1 nn utfk = 2 1 1 1, , 2 2n n k f t h u hk⎛ ⎞= + +⎜ ⎟⎝ ⎠ 3 2 1 1, , 2 2n n k f t h u hk⎛ ⎞= + +⎜ ⎟⎝ ⎠ ).,( 34 hkuhtfk nn ++= （5.22） Butcher表分别为： 常用的四阶Runge-Kutta方法： 100 2 10 2 1 6 1 3 1 3 1 6 1 1 2 1 2 1 0 经典四阶Runge-Kutta方法： ),33( 8 43211 kkkkhuu nn ++++=+ ),,(1 nn utfk = 2 1 1 1, , 3 3n n k f t h u hk⎛ ⎞= + +⎜ ⎟⎝ ⎠ 3 1 2 2 1, , 3 3n n k f t h u hk hk⎛ ⎞= + − +⎜ ⎟⎝ ⎠ ).,( 34 hkuhtfk nn ++= （5.23） Butcher表分别为： 以上讨论的是m级Runge-Kutta法在m=1,2,3,4时，可分别 得到最高阶级一、二、三、四阶，但是，通常m级Runge-Kutta 方法最高阶不一定是m阶。若设p(m)是m级Runge-Kutta方法可 达到的最高阶，可证： (5) 4, (6) 5, (7) 6, (8) 6, (9) 7p p p p p= = = = = 。 1 3 2 3 1 0 100 1 3 1 3 1 − 8 1 8 3 8 3 8 1 1 2 21 , 1 n n n n n t uu u h t+ ⎛ ⎞= + −⎜ ⎟+⎝ ⎠ 改进的Euler法计算公式为： Runge-Kutta法计算公式为： (5.22)求解初值问题: 例1 分别用Euler法，改进的Euler法和Runge-Kutta法 2 21 , 1 tuu t ′ = − + (0) 0u =0 2 ,t≤ ≤ 解：Euler法计算公式为： 1 1 2 1 ( ) , 2n n u u h k k+ = + + 1 221 ,1 n n n t uk t = − + 1 2 2 2( )( )1 , 1 ( ) n n n t h u hkk t h + += − + + ),22( 6 43211 kkkkhuu nn ++++=+ 1 2 21 , 1 n n n t uk t = − + ,) 2 1(1 ) 2 1)( 2 1(2 1 2 1 2 ht hkuht k n nn ++ ++ −= 2 3 2 1 12 2 21 , 1 ( ) n n n t h u hk k t h ⎛ ⎞⎛ ⎞+ +⎜ ⎟⎜ ⎟⎝ ⎠⎝ ⎠= − + + .)(1 ))((21 2 3 4 ht hkuhtk n nn ++ ++−= 三个方法计算结果比较表 Euler法 改进Euler法 Runge-Kutta法 数值解 un 误差 数值解 un 误差 数值解 un 误差 0 0 0 0 0 0 0 0 0 1 0.5 0.433333 0.500000 0.066667 0.400000 0.033333 0.433218 0.000115 2 1.0 0.666667 0.800000 0.133333 0.635000 0.031667 0.666312 0.000355 3 1.5 0.807692 0.900000 0.092308 0.787596 0.020096 0.807423 0.000269 4 2.0 0.933353 0.985615 0.051282 0.921025 0.012308 0.933156 0.000171 n tn 精确解 u(tn) 作比较 ，计算结果见下表： 取步长h=0.5,tn=0.5n,n=0,1,2,3。并与精确解： )1(3 )3()( 2 2 t tttu + += 例如，对最简单的Euler法 L2,1,0,1 =+=+ nhfuu nnn 用其求解模型方程得到 L2,1,0,)1(1 =+=+=+ nuhuhuu nnnn μμ 当un有舍入误差时，其近似解为 nu~ ，从而有 nn uhu ~)1(~ 1 μ+=+ 取 n n ne u u= − % ，得到误差传播方程 1 (1 ) ,n ne h eμ+ = + 单步法的绝对稳定性 记 hh μ= ，只要 11 <+h 差都不会恶性发展，此时方法绝对稳定。 ，则显式Euler方法的解和误 从 1 1,h+ < 可得 2 0h− < < 。即 20 h μ< < − 时， (-1,0)为圆心，1为半径的单位圆。 又由于实数m＜0， 绝对稳定。 若m为复数，在 hh μ= 的复平面上，则 11 <+h 表示为以 2− 0 1− 1 ( )0,1− 绝 对 稳 定 区 域 绝对稳定区间 考虑隐式单步方法中最简单的隐式Euler法 1 1 1( , ) 0, 1,n n n nu u h f t u n+ + += + = L 其特征方程为： ( ) ( )hρ λ σ λ− = (1 ) 1 0h λ− − = 得 1 1 ,1 hλ = − 当 11 >− h 时， 1 1 ,λ < 故 11 >−h 就是 隐式Euler法的绝对稳定区域。 当m＜0为实数时，绝对稳定区间为 (-∞,0)。 h 平面上以（1.0）为圆心的单位圆外区域。它是 当Rem＜0时，它位于 平面上y轴左侧区域。h 1 1h− ≤ 1 1h− > 0Re <μ [ ]2,0∉h ( )0,∞− 隐式Euler法的绝对稳定区域 又如，梯形法 L,1,0)(2 1 11 =++= ++ nffhuu nnnn 其特征方程为： ( ) ( )hρ λ σ λ− = 其根 1( )hλ = 当Rel＜0时， 1 2 1 , 1 2 h h + < − 故梯形公式 h 平面的左半平面。绝对稳定区间为(-∞，0)。的绝对稳定域是 1 2 , 1 2 h h + − 1 1 0 2 2 h hλ⎛ ⎞ ⎛ ⎞− − + =⎜ ⎟ ⎜ ⎟⎝ ⎠ ⎝ ⎠ ( ) 0,1,n m nu P h u nμ= + = L 下面考察Runge-Kutta法的绝对稳定性。 根据定义，在 m=p阶Runge-Kutta法中取 f=mu (m常数)，则 nuk μ=1 2 21(1 ) nk b h uμ μ= + 2 3 3 1 ( )n j j j k u h b kμ = = + ∑ LL nmm uhPk )(1 μμ −= （其中Pi(l)是 i 次多项式）,从而有: 1 1 1 ( ( )) m n n i i n i u u h c P h uμ μ+ − = = + ∑ 1( ) nP h uμ μ= 2 ( ) nP h uμ μ=31 32 1(1 ( )) nb h b hP h uμ μ μ μ= + + 注意， uu μ=′ 的解 uetu μ=)( 且 pjutu jj ,,1,0)()( L== μ 则 ( )1)(21 )(!)(!2)()()( ++ +++′′+′+= pnp p nnnn hOtup htuhtuhtutu L ( )12 )() ! )( !2 )(1( ++++++= pn p hOtu p hhh μμμ L 若为p阶方法，则应有 ( ) ( )( ) )()(1 ! )( 1 0 + = +⎥⎦ ⎤⎢⎣ ⎡ +−= ∑ pnp k m k hOtuhP k hhR μμ 从而 ⎟⎟⎠ ⎞ ⎜⎜⎝ ⎛=+ ∑ = p k k m k hhP 0 ! )()(1 μμ 记 hh μ= ，则可将上式写成 L,1,0)( !0 1 ==⎟⎟⎠ ⎞ ⎜⎜⎝ ⎛= ∑ = + nuhuk hu nn m k k n λ 进而误差传播方程为： nn h ελε )(1 =+ 其中 . ! 1 !2 11)( 2 mh m hhh ++++= Lλ 注意，当 4,3,2,1=m 时，解不等式 )(hλ <1就可得显式 4,3,2,1 Runge-Kutta法公式绝对稳定域。 当m＜0为实数，则得各阶 （ m=1,2,3,4 ）的绝对稳定区间(见下表)。 =m 级 绝对稳定区间 一级 （-2，0） 二级 （-2，0） 三级 （-2.51，0） 四级 （-2.78，0） 32 6 1 !2 11 hhh +++ 2 !2 11 hh ++ h+1 )(hλ 432 24 1 6 1 2 11 hhhh ++++ Runge-Kutta法（ m=1,2,3,4 ）的绝对稳定区间表 使用Runge-Kutta法时，应取h属于绝对稳定区域，否则有 可能产生很大误差。 例如用四阶Runge-Kutta法（5.22）求解 ,20uu −=′ ( )0 1u = 取步长h=0.1及0.2。当步长为0.1时， 属于绝对稳定区域，误差h 随n下降。 不属于绝对稳定区域，误差增长很当步长为0.2时，h 快。 t h=0.1的误差 h=0.2的误差 0.0 0.0 0.0 0.2 -0.092 795 4.98 0.4 -0.012 010 2.50 0.6 -0.001 366 125.0 0.8 -0.000 152 625.0 1.0 -0.000 017 3125.0 数值计算结果误差表 THE END