第7章-2 收斂性穩(wěn)定性R-K方法

1、常微分方程的數(shù)值解法的常微分方程的數(shù)值解法的 收斂性、穩(wěn)定性收斂性、穩(wěn)定性第第7 7章章-2-2 以上我們討論了求解問題（以上我們討論了求解問題（7-1）, ,（7-2）的單步）的單步法法和多步法。和多步法。應(yīng)關(guān)注三個問題：應(yīng)關(guān)注三個問題：、數(shù)值方法的局部截斷誤差和階、數(shù)值方法的局部截斷誤差和階二、在離散點二、在離散點tn處的數(shù)值解處的數(shù)值解un是否收斂到精確解是否收斂到精確解u( (tn) ) 三、數(shù)值方法的穩(wěn)定性三、數(shù)值方法的穩(wěn)定性具體說，具體說，對于上述兩類方法求近似解（數(shù)值解）還對于上述兩類方法求近似解（數(shù)值解）還誤差估計誤差估計、收斂性和穩(wěn)定性收斂性和穩(wěn)定性。 對于第一個問題前面我們

2、已經(jīng)討論過，而關(guān)于數(shù)值對于第一個問題前面我們已經(jīng)討論過，而關(guān)于數(shù)值方法收斂性問題我們在這里不詳細討論，只給出一些基方法收斂性問題我們在這里不詳細討論，只給出一些基本結(jié)論性的結(jié)果，即：本結(jié)論性的結(jié)果，即：對單步法，當方法的階對單步法，當方法的階p11時，有整體誤差時，有整體誤差)()(pnnnhOutuE故有故有 0lim0nhE，因此方法是收斂的。，因此方法是收斂的。 對于多步法，若方法是對于多步法，若方法是k 步步p 階法，那么（階法，那么（7-24）是）是 一個一個k階差分方程，引入多步法（階差分方程，引入多步法（7-24）的）的第一特征第一特征多多項項式式和和第二特征多項式第二特征多項式

3、： 0( ),kjjj kjjj0)(定義定義7.17.1 若（若（7-24）的第一特征多項式）的第一特征多項式( () )的所有的所有 根在單位圓內(nèi)或圓上（根在單位圓內(nèi)或圓上（11），且位于單位圓周上），且位于單位圓周上的根都是單根，稱多步法（的根都是單根，稱多步法（7-24）滿足）滿足根條件根條件。第二特征多項式第二特征多項式第一特征多項式第一特征多項式定理定理7.27.2 若線性多步法若線性多步法( (7-24) )的階的階p11，且滿足，且滿足根條件，則方法是收斂的。根條件，則方法是收斂的。 對于常用的數(shù)值方法都是滿足收斂性條件的對于常用的數(shù)值方法都是滿足收斂性條件的。 下面我們著重討

4、論第三個問題，即下面我們著重討論第三個問題，即數(shù)值方法的穩(wěn)數(shù)值方法的穩(wěn)是有誤差的，且這些誤差將在計算中傳遞下去。是有誤差的，且這些誤差將在計算中傳遞下去。定性定性問題。問題。110,kuuu誤差積累無限增長，則會歪曲真解，這樣的算法是不誤差積累無限增長，則會歪曲真解，這樣的算法是不如果如果能用的。能用的。用多步法計算時，各種因素如初值用多步法計算時，各種因素如初值精確解為精確解為 2 2( ) (1)u tt ?？紤]二步考慮二步三階顯式法三階顯式法： 21145(42)nnnnnuuuhff例如例如 初值問題初值問題124,utu 02;t (0)1u取步長取步長h=0.1=0.1，初值初值u

5、0=1=1，附加值：附加值： 2 21(1) (0.1)uhh。數(shù)值結(jié)果表數(shù)值結(jié)果表在開始幾步數(shù)值解與精確解符合，但在再往后算，數(shù)值在開始幾步數(shù)值解與精確解符合，但在再往后算，數(shù)值解的誤差則急劇增長，完全歪曲了真解解的誤差則急劇增長，完全歪曲了真解. .通常人們都是通過模型方程來討論方法的通常人們都是通過模型方程來討論方法的數(shù)值穩(wěn)定性數(shù)值穩(wěn)定性。uu( (7-32) )而一般形式的一階微分方程總能化成（而一般形式的一階微分方程總能化成（7-32）的形式。）的形式。因為實際計算時，。因為實際計算時，h是固定的。是固定的。 當某一步當某一步un有舍入誤差時，有舍入誤差時， 若以后的計算中不會逐步擴

6、大，稱這種穩(wěn)定性為若以后的計算中不會逐步擴大，稱這種穩(wěn)定性為絕對穩(wěn)定性絕對穩(wěn)定性。此后，若不做特殊說明，都是指絕對穩(wěn)定性此后，若不做特殊說明，都是指絕對穩(wěn)定性 。模型方程為：模型方程為：本書中數(shù)值方法的穩(wěn)定性也是如此。前提是求解好條件問題，本書中數(shù)值方法的穩(wěn)定性也是如此。前提是求解好條件問題，其中其中Re( () )0 0。另外，我們也不考慮另外，我們也不考慮h00時方法的漸近穩(wěn)定性時方法的漸近穩(wěn)定性 例如例如，對最簡單的，對最簡單的Euler法法2, 1, 0,1nhfuunnn( (7-33) ) 用其求解模型方程（用其求解模型方程（7-32）得到）得到1nnnuuh u 當當un有舍入誤

7、差時，其近似解為有舍入誤差時，其近似解為 nu，從而有，從而有 nnuhu)1(1取取 nnnuu，得到誤差傳播方程，得到誤差傳播方程 ,)1 (1nnh(1),0,1, 2nh un記記 hh，只要，只要 11h都不會惡性發(fā)展，此時方法絕對穩(wěn)定。都不會惡性發(fā)展，此時方法絕對穩(wěn)定。，則顯式，則顯式Euler方法的解和誤差方法的解和誤差從從 11,h可得可得 20h 。即即20 h 時，時，(-1,0)(-1,0)為圓心，為圓心，1 1為半徑的單位圓。為半徑的單位圓。 又由于實數(shù)又由于實數(shù)0 0，（7-33）絕對穩(wěn)定，絕對穩(wěn)定，若若為復(fù)數(shù)，在為復(fù)數(shù)，在 hh的復(fù)平面上，則的復(fù)平面上，則 11h表

8、示為以表示為以20110, 1絕絕對對穩(wěn)穩(wěn)定定區(qū)區(qū)域域絕對穩(wěn)定區(qū)間絕對穩(wěn)定區(qū)間 定義定義7.27.2 一個數(shù)值方法用于求解模型問題（一個數(shù)值方法用于求解模型問題（7-32），若在），若在 平面中的某一區(qū)域平面中的某一區(qū)域D中方法都是絕對穩(wěn)定的，而在區(qū)域中方法都是絕對穩(wěn)定的，而在區(qū)域D外，方法外，方法是不穩(wěn)定的是不穩(wěn)定的，則稱則稱D是方法的是方法的絕對穩(wěn)定區(qū)域；絕對穩(wěn)定區(qū)域；絕對穩(wěn)定區(qū)間絕對穩(wěn)定區(qū)間。它與實軸的交稱為它與實軸的交稱為例如，顯式例如，顯式Euler方法的方法的絕對穩(wěn)定區(qū)域、區(qū)間。如圖絕對穩(wěn)定區(qū)域、區(qū)間。如圖現(xiàn)在考察多步法現(xiàn)在考察多步法( (7-24) )，將它用于解模型方程，將它用

9、于解模型方程( (7-32) )得到得到k階線性差分方程階線性差分方程0kjnjju(7-34) 若取若取 hh，則記（，則記（7-34）的特征方程為）的特征方程為0)()(h(7-35) 其中其中 kjjj0)(kjjj0)(0kjnjjhu由由k階線性差分方程的性質(zhì)我們可以得到如下結(jié)論，階線性差分方程的性質(zhì)我們可以得到如下結(jié)論，區(qū)域：區(qū)域： ()1,1,2,jhhjkD例如，對于例如，對于k=1=1時，考慮隱式方法中最簡單的后退時，考慮隱式方法中最簡單的后退Euler法法111(,)0, 1,nnnnuuh f tun方程（方程（7-35）的根都在單位圓內(nèi)（）的根都在單位圓內(nèi)（1) 1)

10、，則線性多步法，則線性多步法（7-4）關(guān)于）關(guān)于 hh絕對穩(wěn)定絕對穩(wěn)定, ,其絕對穩(wěn)定域是復(fù)平面其絕對穩(wěn)定域是復(fù)平面 h上的上的其特征方程為：其特征方程為：( )( )h (1)10h 若特征若特征得得 11,1h當當 11 h時，時， 11 ,故故 11h就是就是隱式隱式Euler法的絕對穩(wěn)定區(qū)域。法的絕對穩(wěn)定區(qū)域。當當0 0為實數(shù)時，絕對穩(wěn)定區(qū)間為為實數(shù)時，絕對穩(wěn)定區(qū)間為 (-,0)(-,0)。h平面上以（平面上以（1.01.0）為圓心的單位圓外區(qū)域。）為圓心的單位圓外區(qū)域。它是它是11h11h0Re2,0h0,當當Re 0 0時，它位于時，它位于 平面上平面上y軸左側(cè)區(qū)域軸左側(cè)區(qū)域。h又

11、如，梯形法又如，梯形法, 1 , 0)(2111nffhuunnnn其特征方程為：其特征方程為：( )( )h 其根其根1( )h當當Re0 0時，時，121,12hh故梯形公式故梯形公式h平面的左半平面。絕對穩(wěn)定區(qū)間為平面的左半平面。絕對穩(wěn)定區(qū)間為(-(-，0 0) )。的絕對的絕對穩(wěn)定域是穩(wěn)定域是 這樣檢驗絕對穩(wěn)定性歸結(jié)為檢驗特征方程這樣檢驗絕對穩(wěn)定性歸結(jié)為檢驗特征方程( (7-35) )的根是否在單位的根是否在單位圓內(nèi)（圓內(nèi)（1)1)。對此有對此有很多判別法，如很多判別法，如Schur準則、軌跡法。準則、軌跡法。12,12hh11022hhk=1=14 4的隱式的隱式Adams類方法的絕

12、對穩(wěn)定區(qū)間（類方法的絕對穩(wěn)定區(qū)間（0 0為實數(shù)）。為實數(shù)）。實系數(shù)二次方程實系數(shù)二次方程2 2- -b - -c=0=0的根在單位圓內(nèi)的充要條件為的根在單位圓內(nèi)的充要條件為: : 21cb這里我們給出一種簡單的、常用的判別法：這里我們給出一種簡單的、常用的判別法：例例 證明求解一階常微分方程初值問題：證明求解一階常微分方程初值問題：),( utfu 0)0(uu的差分格式的差分格式)85 (121212nnnnnfffhuu收斂收斂并求其并求其局部截斷誤差主項局部截斷誤差主項、絕對穩(wěn)定區(qū)間絕對穩(wěn)定區(qū)間。解解：由差分格式可知，：由差分格式可知，2( ), ()(1)0 , 則則其特征值其特征值滿

13、足根條件滿足根條件。令令得得1=0,=0,2=1=1。2581( ),121212 故此為故此為隱式二步三階隱式二步三階法，其局部截斷誤差主項為：法，其局部截斷誤差主項為：4(4)1()24nh ut。01 10C 1C 2C 3C 4C 1852 1121212 013( 14)22 032118512262 1212 776604311851224!3! 12120注意，注意，0120,1,1,從而從而012185,121212由定理由定理7.27.2 可知，此方法收斂?？芍朔椒ㄊ諗?。( )( )h 12 812 5hh而而 1212 5hh自然成立。自然成立

14、。hhhh512412512812得得 12 4h 12 812 4hh即有即有 332,hh 可得其可得其絕對穩(wěn)定區(qū)間絕對穩(wěn)定區(qū)間： 60h 。又其特征方程為又其特征方程為21280125125hhhh而使得而使得 1 1的充要條件為：的充要條件為：124125hh12125hh現(xiàn)在再由現(xiàn)在再由128 h124 h2581110121212hhh進一步進一步而而 自然成立。自然成立。 顯式顯式Runge-KuttaRunge-Kutta法法第第7 7章章-3-37.1.47.1.4 四階顯式四階顯式Runge-Kutta法法 通過觀察我們發(fā)現(xiàn)顯式通過觀察我們發(fā)現(xiàn)顯式Euler法和隱法和隱Eu

15、ler法各用到了法各用到了u( (t) )在在 t,t+h 上的一個一階導(dǎo)數(shù)值，它們都是一階方法。上的一個一階導(dǎo)數(shù)值，它們都是一階方法。改進的改進的Euler法用到了法用到了u( (t) )在在 t,t+h 上的兩個一階導(dǎo)數(shù)值，它上的兩個一階導(dǎo)數(shù)值，它 梯形法和梯形法和們都是二階方法。們都是二階方法。 我們要研究的我們要研究的Runge-Kutta方法是一種高階單步法，它使用方法是一種高階單步法，它使用u( (t) )在在 t,t+h 上的斜率上的斜率f 在一些點的值非線性表示在一些點的值非線性表示 使得其局部截斷誤差的階和使得其局部截斷誤差的階和Taylor展開法相等。展開法相等。( , (

16、 ), )t u t h Euler是最簡單的單步法。單步法不需要附加初值，所需的是最簡單的單步法。單步法不需要附加初值，所需的存儲量小，改變步長靈活，但線性單步法的階最高為存儲量小，改變步長靈活，但線性單步法的階最高為2 2，Taylor展開法，用在同一點展開法，用在同一點( (tn,un) )的高階導(dǎo)數(shù)表示的高階導(dǎo)數(shù)表示 ，這不便，這不便于計算。于計算。( , ( ), )t u t h先引進若干記號，首先先引進若干記號，首先 t, t+h 取上的取上的m個點個點: :123,mtttttth令令, 2,1mihathattiii 11,2,iijijbaimmiiicc11, 02131

17、32121mmmmbbbbbb Runge-Kutta矩陣矩陣B為嚴格下三角矩陣：為嚴格下三角矩陣：滿足滿足顯式顯式 Runge-Kutta Runge-Kutta 公式公式假設(shè)三組系數(shù)已給定，則求解（假設(shè)三組系數(shù)已給定，則求解（7-1），（），（7-2）的一般）的一般1( , ),nnnnuuht u h （7-12）其中其中1( , )mnniiituhc k( (7-13) )( (7-14) )顯式顯式Runge-Kutta法的計算過程如下：法的計算過程如下：1( ,),nnkf t u2221 1(,) ,nnkf tha uhb k3331 1322(,(),nnkf tha uh

18、 b kb k31323,bba11(,),mmnmnmjjjkf thauhb k 212,ba11mmjmjba0,1,n 現(xiàn)在推導(dǎo)一些常用的計算方案，特別地，給出現(xiàn)在推導(dǎo)一些常用的計算方案，特別地，給出 m=3 =3 顯式顯式首先將首先將u( (t+h) )在在t處展開到處展開到h的三次冪的三次冪, ,即：即：3( )41()( )( )()!lllhu thu tutO hl( (7-15) )其中其中( (7-16) ) Runge-Kutta法的推導(dǎo)。法的推導(dǎo)。( )( , ( ), ),u tht u t h2311( , , )()()26ut u hfhfhfffO htuf

19、fff22tttuuuffffff其次，由二元函數(shù)其次，由二元函數(shù)f(t,u(t)在在( (t, ,u) )點處的點處的Taylor展開式可得：展開式可得： ,)(,(1ftutfk),(1222khauhatfk21()tufha fk f)(2132222hOfahfhaf)()21(323322233hOfaffbahfhafku于是，將于是，將k1 1, ,k2 2, ,k3 3代入（代入（7-13）中，即）中，即31( , , )iiit u hc k1c f22332 32312ucfa f ha b f fa fh1 12233c kc kc k222321112()2tttuu

20、uh afk fk fO h222312cfha fa f h3()O h1232 23 3cccfh a ca cf31( , , )iiit u hck22232 32 3223312()()2uha b c ffa ca cfO h( (7-17) ) 由由( (7-16) )已得已得(3)211( , , )( )( )( )26t u hu tu th uth其中其中 ,tuffff22tttuuufffff f。2311()26ufh fhfffO h合并合并f( (t, ,u, ,h) )展開式中的各階展開式中的各階hl( (l= =0,1,2 2) )的系數(shù)，得的系數(shù)，得比較比

21、較 ( , , )t u h和和 ( , , )t u h的同次冪系數(shù)，的同次冪系數(shù)，231111()2233ufh fhfffO h可得可得（一）（一）m=1 =1 此時此時 c2= =c3= =0, , f( (t, ,u, ,h)=)=c1 f , , 比較比較h的零次冪，知的零次冪，知( , , ),t u hf方法（方法（7-21）為一級一階）為一級一階Runge-Kutta法，實際上為法，實際上為Euler法。法。 （二）（二）m=2,=2,此時此時 c3=0，則則223122 2221( , , )()2t u hccfha c fh a c fO h與與 ( , , )t u

22、h比較比較1, ,h的系數(shù)，則的系數(shù)，則 2112221cacc它有無窮多組解，從而有無窮多個二級二階方法。它有無窮多組解，從而有無窮多個二級二階方法。2311()26ufh fhfffO h( , , )t u h,21, 1, 0221acc（1 1）稱為稱為中點法中點法。此時此時12,nnuuhk1( ,),nnkf t u2111,22nnkfth uhk(7-18) 三個常見的方法是：三個常見的方法是：1221,1,2cca（2 2）稱為稱為改進的改進的Euler法法。此時此時112,2nnhuukk1( ,),nnkf t u21(,)nnkf th uhk(7-19) 12213

23、2,443cca（3 3）此時此時1123,4nnhuukk1( ,),nnkf t u2122,33nnkfth uhk（三）（三）m=3 3 比較（比較（7-16）和（）和（7-17），令），令 f,f,h, ,h2 2 的系數(shù)的系數(shù)相等，并注意相等，并注意 ff,的任意性，得的任意性，得 21, 13322321cacaccc.61,313322323222cbacaca四個方程不能完全確定六個系數(shù)，因此這是含兩個參數(shù)的四個方程不能完全確定六個系數(shù)，因此這是含兩個參數(shù)的三級三級三階三階方法類。方法類。常見方案有：常見方案有：12313,0,44ccc2332122,333aab。 Heu

24、n三階方法三階方法。 此時取此時取1133,4nnhuukk1( ,),nnkf t u2111,33nnkfth uhk(7-20) 3222,33nnkfth uhk（2 2）Kutta三階方法三階方法，123121,636ccc23321,1,22aab 。),4(613211kkkhuunn),(1nnutfk ),21,21(12hkuhtfknn).2,(213hkhkuhtfknn(7-21) (7-21) 此時此時（四）（四）m=4=4將（將（7-167-16）和（）和（7-177-17）展開到）展開到h3, ,比較比較)3 , 2 , 1 , 0( ihi的系數(shù)，則含的系數(shù)，

25、則含1313個待定系數(shù)的個待定系數(shù)的1111個方程，由此得到含兩個參數(shù)的個方程，由此得到含兩個參數(shù)的四級四階四級四階Runge-Kutta方法方法類，類，其中最常用的有以下兩個方法：其中最常用的有以下兩個方法： 經(jīng)典四階經(jīng)典四階Runge-Kutta方法方法：),22(643211kkkkhuunn),(1nnutfk 2111,22nnkf th uhk3211,22nnkf th uhk).,(34hkuhtfknn（7-22）),33(843211kkkkhuunn),(1nnutfk 2111,33nnkf th uhk31221,33nnkf th uhkhk).,(34hkuhtf

26、knnButcher表分別為：表分別為：10021021613131611212101001313181838381132310以上討論的是以上討論的是m級級Runge-Kutta法在法在m=1,2,3,4=1,2,3,4時，可分別時，可分別得到最高階級一、二、三、四階，但是，通常得到最高階級一、二、三、四階，但是，通常m級級Runge-Kutta Runge-Kutta 方法最高階不一定是方法最高階不一定是m階。若設(shè)階。若設(shè)p( (m) )是是m級級Runge-Kutta方法可方法可達到的最高階，可證達到的最高階，可證：(5)4,(6)5,(7)6,(8)6,(9)7ppppp 。1221,

27、1nnnnnt uuuht改進的改進的Euler法計算公式為：法計算公式為：經(jīng)典經(jīng)典Runge-Kutta法計算公式為：法計算公式為：例例1 分別用分別用EulerEuler法，改進的法，改進的EulerEuler法法(7-27)(7-27)和經(jīng)典和經(jīng)典Runge-Kutta法法(7-30)(7-30)求解初值問題求解初值問題: :221,1tuut (0)0u02,t 解解：Euler法計算公式為：法計算公式為： 1121(),2nnuuh kk1221,1nnnt ukt 1222()()1,1 ()nnnth uhkkth ),22(643211kkkkhuunn1221,1nnnt ukt ,)21(1)21)(21( 21212hthkuhtknnn232112221,1 ()nnnthuhkkth .)(1)(21234hthkuhtknnn三個方法計算結(jié)果比較表三個方法計算結(jié)果比較表作比較作比較 ，計算結(jié)果見下表：，計算結(jié)果見下表：取步長取步長h=0.5,=0.5,tn=0.5=0.5n, ,n