第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)注三個(gè)問題：應(yīng)關(guān)注三個(gè)問題：、數(shù)值方法的局部截?cái)嗾`差和階、數(shù)值方法的局部截?cái)嗾`差和階二、在離散點(diǎn)二、在離散點(diǎn)tn處的數(shù)值解處的數(shù)值解un是否收斂到精確解是否收斂到精確解u( (tn) ) 三、數(shù)值方法的穩(wěn)定性三、數(shù)值方法的穩(wěn)定性具體說，具體說，對(duì)于上述兩類方法求近似解（數(shù)值解）還對(duì)于上述兩類方法求近似解（數(shù)值解）還誤差估計(jì)誤差估計(jì)、收斂性和穩(wěn)定性收斂性和穩(wěn)定性。 對(duì)于第一個(gè)問題前面我們

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

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

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

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

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

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

8、示為以表示為以20110, 1絕絕對(duì)對(duì)穩(wěn)穩(wěn)定定區(qū)區(qū)域域絕對(duì)穩(wěn)定區(qū)間絕對(duì)穩(wěn)定區(qū)間 定義定義7.27.2 一個(gè)數(shù)值方法用于求解模型問題（一個(gè)數(shù)值方法用于求解模型問題（7-32），若在），若在 平面中的某一區(qū)域平面中的某一區(qū)域D中方法都是絕對(duì)穩(wěn)定的，而在區(qū)域中方法都是絕對(duì)穩(wěn)定的，而在區(qū)域D外，方法外，方法是不穩(wěn)定的是不穩(wěn)定的，則稱則稱D是方法的是方法的絕對(duì)穩(wěn)定區(qū)域；絕對(duì)穩(wěn)定區(qū)域；絕對(duì)穩(wěn)定區(qū)間絕對(duì)穩(wěn)定區(qū)間。它與實(shí)軸的交稱為它與實(shí)軸的交稱為例如，顯式例如，顯式Euler方法的方法的絕對(duì)穩(wěn)定區(qū)域、區(qū)間。如圖絕對(duì)穩(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例如，對(duì)于例如，對(duì)于k=1=1時(shí)，考慮隱式方法中最簡(jiǎn)單的后退時(shí)，考慮隱式方法中最簡(jiǎn)單的后退Euler法法111(,)0, 1,nnnnuuh f tun方程（方程（7-35）的根都在單位圓內(nèi)（）的根都在單位圓內(nèi)（1) 1)

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

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

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

13、足根條件滿足根條件。令令得得1=0,=0,2=1=1。2581( ),121212 故此為故此為隱式二步三階隱式二步三階法，其局部截?cái)嗾`差主項(xiàng)為：法，其局部截?cái)嗾`差主項(xiàng)為：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 可得其可得其絕對(duì)穩(wěn)定區(qū)間絕對(duì)穩(wěn)定區(qū)間： 60h 。又其特征方程為又其特征方程為21280125125hhhh而使得而使得 1 1的充要條件為：的充要條件為：124125hh12125hh現(xiàn)在再由現(xiàn)在再由128 h124 h2581110121212hhh進(jìn)一步進(jìn)一步而而 自然成立。自然成立。 顯式顯式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 上的一個(gè)一階導(dǎo)數(shù)值，它們都是一階方法。上的一個(gè)一階導(dǎo)數(shù)值，它們都是一階方法。改進(jìn)的改進(jìn)的Euler法用到了法用到了u( (t) )在在 t,t+h 上的兩個(gè)一階導(dǎo)數(shù)值，它上的兩個(gè)一階導(dǎo)數(shù)值，它 梯形法和梯形法和們都是二階方法。們都是二階方法。 我們要研究的我們要研究的Runge-Kutta方法是一種高階單步法，它使用方法是一種高階單步法，它使用u( (t) )在在 t,t+h 上的斜率上的斜率f 在一些點(diǎn)的值非線性表示在一些點(diǎn)的值非線性表示 使得其局部截?cái)嗾`差的階和使得其局部截?cái)嗾`差的階和Taylor展開法相等。展開法相等。( , (

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

17、32121mmmmbbbbbb Runge-Kutta矩陣矩陣B為嚴(yán)格下三角矩陣：為嚴(yán)格下三角矩陣：滿足滿足顯式顯式 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法的計(jì)算過程如下：法的計(jì)算過程如下：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)一些常用的計(jì)算方案，特別地，給出現(xiàn)在推導(dǎo)一些常用的計(jì)算方案，特別地，給出 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) )點(diǎn)處的點(diǎn)處的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 此時(shí)此時(shí) c2= =c3= =0, , f( (t, ,u, ,h)=)=c1 f , , 比較比較h的零次冪，知的零次冪，知( , , ),t u hf方法（方法（7-21）為一級(jí)一階）為一級(jí)一階Runge-Kutta法，實(shí)際上為法，實(shí)際上為Euler法。法。 （二）（二）m=2,=2,此時(shí)此時(shí) c3=0，則則223122 2221( , , )()2t u hccfha c fh a c fO h與與 ( , , )t u

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

23、2,443cca（3 3）此時(shí)此時(shí)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四個(gè)方程不能完全確定六個(gè)系數(shù)，因此這是含兩個(gè)參數(shù)的四個(gè)方程不能完全確定六個(gè)系數(shù)，因此這是含兩個(gè)參數(shù)的三級(jí)三級(jí)三階三階方法類。方法類。常見方案有：常見方案有：12313,0,44ccc2332122,333aab。 Heu

24、n三階方法三階方法。 此時(shí)取此時(shí)取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) 此時(shí)此時(shí)（四）（四）m=4=4將（將（7-167-16）和（）和（7-177-17）展開到）展開到h3, ,比較比較)3 , 2 , 1 , 0( ihi的系數(shù)，則含的系數(shù)，

25、則含1313個(gè)待定系數(shù)的個(gè)待定系數(shù)的1111個(gè)方程，由此得到含兩個(gè)參數(shù)的個(gè)方程，由此得到含兩個(gè)參數(shù)的四級(jí)四階四級(jí)四階Runge-Kutta方法方法類，類，其中最常用的有以下兩個(gè)方法：其中最常用的有以下兩個(gè)方法： 經(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級(jí)級(jí)Runge-Kutta法在法在m=1,2,3,4=1,2,3,4時(shí)，可分別時(shí)，可分別得到最高階級(jí)一、二、三、四階，但是，通常得到最高階級(jí)一、二、三、四階，但是，通常m級(jí)級(jí)Runge-Kutta Runge-Kutta 方法最高階不一定是方法最高階不一定是m階。若設(shè)階。若設(shè)p( (m) )是是m級(jí)級(jí)Runge-Kutta方法可方法可達(dá)到的最高階，可證達(dá)到的最高階，可證：(5)4,(6)5,(7)6,(8)6,(9)7ppppp 。1221,

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