1、常微分方程的數(shù)值解法常微分方程的數(shù)值解法第第7 7章章0( )( , ( )dtau tufuuuLutfutf),(),(考慮常微分方程的初值問題考慮常微分方程的初值問題 (7-1) (7-2)則（則（7-1）的解存在且唯一。）的解存在且唯一。 或與其等價的積分方程或與其等價的積分方程 若若f( (t,u) )關(guān)于關(guān)于u u滿足滿足Lipschitz條件，即存在常數(shù)條件，即存在常數(shù)L，對任意，對任意ta, ,b,均有均有( , ),uf t uatb 0( )u au 它是一種離散化方法，利用這種方法，可以在一系列事它是一種離散化方法，利用這種方法，可以在一系列事先先btttaN21)(,)

2、,(),(21Ntututu求出其求出其上的上的未知函數(shù)未知函數(shù)u( (t) )之值之值的近似值的近似值: :Nuuu,21而而u1, u2, , uN通常稱為通常稱為初值問題的數(shù)值解初值問題的數(shù)值解。 首先我們利用數(shù)值積分公式建立求解（首先我們利用數(shù)值積分公式建立求解（7-17-1）或（）或（7-27-2）的數(shù)）的數(shù)值方法。值方法。 什么是什么是數(shù)值解法數(shù)值解法？通常取成等距，即通常取成等距，即ti=t0+ih, i=1, 2, , N,其中其中h0稱為步長。稱為步長。取定的取定的 a, ,b 中的離散點（稱為節(jié)點）中的離散點（稱為節(jié)點）基于數(shù)值積分的解法基于數(shù)值積分的解法 ,ntanh,N

3、abhNn, 2, 1, 0由（由（7-2）式）式將節(jié)點取為將節(jié)點取為（7-3） 如果如果u(tn)的近似值的近似值un已經(jīng)求出，則通過計算已經(jīng)求出，則通過計算(7-3)(7-3)右端右端項的數(shù)值積分可求出項的數(shù)值積分可求出u( (tn+1) )的近似值的近似值un+1+111()( )( , ( )nntnntu tu tf t u t dt0( )u au1)(,()()(1nnttnndttutftutu),(1nnnnutfhuu1, 2, 1, 0Nn令令 上式稱為上式稱為Euler求解公式求解公式，又稱，又稱矩形公式矩形公式。一、一、Euler法法)(,()(nnntutfhtu首

4、先，對首先，對( (7-2) )右端積分項使用左矩形求積公式，右端積分項使用左矩形求積公式，則得則得22( )10001u ttut 的解的解u(t)在在t=0.3=0.3處的數(shù)值解處的數(shù)值解u3 3。 例例：（取步長（取步長h=0.1,小數(shù)點后保留小數(shù)點后保留4 4位）位） 0)0(u解：解： 相應(yīng)的相應(yīng)的EulerEuler公式：公式：221100nnnnuthuu由初值由初值0)0(0 uu，計算得，計算得) 1 . 0(u2) 2 . 0(uu221001 . 0nnnutu1u0 . 01000 . 01 . 00 . 00000. 0202001001 . 0utu0 . 0100

5、1 . 01 . 00 . 020010. 0212111001 . 0utu3) 3 . 0(uu220010. 01002 . 01 . 00010. 00051. 0222221001 . 0utu11()( )( , ( )nntnntu tu tf t u t dt),(111nnnnutfhuu1, 2, 1, 0Nn首先，首先，對對(7-3)右端積分項使用右矩形求積公式，則得右端積分項使用右矩形求積公式，則得令令 上式稱為上式稱為隱式隱式Euler公式公式，又稱，又稱右右矩形公式矩形公式。隱式隱式Euler法法11( )(, ()nnnu th f tu t二、梯形法二、梯形法對

6、對(7-2)右端的積分使用梯形求積分式計算右端的積分使用梯形求積分式計算，則得則得令令111( , )(,)2nnnnnnhuuf t uf tuNn, 2, 1, 0上式稱為上式稱為梯形求解公式梯形求解公式，簡稱，簡稱梯形法梯形法 )(1ntu)(ntu1)(,(nnttdttutf)(,()(,(211nnnntutftutfh)(,()(,(211nnnntutftutfh將將Euler公式與隱式公式與隱式Euler公式做算術(shù)平均，也可得出梯形公式公式做算術(shù)平均，也可得出梯形公式),(),(211nnnnnutfutfhu0121nu123取初值為取初值為 ， nnuu01反復(fù)迭代，即反

7、復(fù)迭代，即 01nu一般的迭代公式表示為：一般的迭代公式表示為：k= 0,1,2,= 0,1,2, ),(),(211knnnnnutfutfhu11knu梯形公式與梯形公式與Euler公式相比要精確的多，但是梯形公式的公式相比要精確的多，但是梯形公式的從而要用如下迭代公式：從而要用如下迭代公式：每步計算要解一個關(guān)于每步計算要解一個關(guān)于un+1+1的非線性方程，的非線性方程，計算量要大一些。計算量要大一些。 01nu 11nu 21nu， ， 若序列若序列 收斂于收斂于 ，當(dāng)，當(dāng) 時，得到：時，得到： 01kknu *1nuk ),(),(2*11*1nnnnnnutfutfhuu則取則取 為

8、第為第 個近似值。個近似值。 1nu *1nu1n如此如此迭代迭代下去下去得到迭代序列得到迭代序列： knu1， knknuu111為終止條件，此時取為終止條件，此時取 作為作為 的近似值的近似值 。 11knu)(1ntu1nu為了避免求解非線性代數(shù)方程，可以用為了避免求解非線性代數(shù)方程，可以用Euler法將它顯化，法將它顯化，0u1uNu1u2u2uNu),(1nnnnuthfuu),(),(2111nnnnnnutfutfhuu00)(utu建立建立預(yù)測預(yù)測校正系統(tǒng)校正系統(tǒng)：此求解公式稱為此求解公式稱為改進的改進的Euler法法，其中，其中 稱為預(yù)測值，稱為預(yù)測值，1nu稱為校正值稱為校

9、正值. . 1nu其求解順序為：其求解順序為：改進的改進的Euler法還可寫成如下形式：法還可寫成如下形式：11( ,)(,)2nnnnnhuuf t uf t如果如果f( (t, ,u( (t)關(guān)于關(guān)于u是線性函數(shù)，則隱式公式可以顯式化。是線性函數(shù)，則隱式公式可以顯式化。 例例 若方程為：若方程為： 5)(uttu隱式隱式Euler公式公式： )5(111nnnnuthuuhthuunnn1115Nn, 2, 1, 0， 梯形公式梯形公式： 102111nnnnnnututhuuhuththunnnn52121111Nn, 2, 1, 0( ,)nnnuh f t u三、三、Milne公式公

10、式若若使用使用 Simpson求積公式，得求積公式，得 dttutfnntt2)(,(令令 ),(,4),(6211222nnnnnnnnutfutfutfhuu進一步可寫成進一步可寫成nnnnnfffhuu12243其中其中 ),(222nnnutff此為二步方法，需要已知此為二步方法，需要已知un和和un+1+1，才能由上式計算出，才能由上式計算出un+2+2的值。的值。),(111nnnutff),(nnnutff )(,(,4)(,(611222nnnnnnnntutftutftutftt二步以上的方法也稱為二步以上的方法也稱為多步法多步法。定義定義 ( )nE h 為求解公式在為求解

11、公式在tn點上的點上的整體截斷誤差整體截斷誤差。衡量求解公式好壞的一個主要標(biāo)準(zhǔn)是求解公式的精度。衡量求解公式好壞的一個主要標(biāo)準(zhǔn)是求解公式的精度。定義定義 假設(shè)假設(shè)ui=(=(u( (ti)，i=0,1,2,=0,1,2, ,n-1，則稱，則稱nnnutuhR)()(為求解公式第為求解公式第n步的步的局部截斷誤差局部截斷誤差。nt1ntt tu11nnu tu nu t1nu( )( )nnnR hu tu nEho( )nnu tu1( )niiR h1nu 如果設(shè)某求解公式的局部截斷誤如果設(shè)某求解公式的局部截斷誤差差：)()(1pnhOhR這樣我們就稱該求解公式具有這樣我們就稱該求解公式具有

12、 p p 階精度階精度。)()(pnhOhE則我們可以證明其整體截斷誤差為：則我們可以證明其整體截斷誤差為：事實上，若事實上，若niinhRhE1)()(niphO11)(niphOh1)(nhOhp)(nabnhOp)()(phO),()(1pihOhR, 2, 1ni則則 求解公式的精度越高，計算解的精確性可能越好。求解公式的精度越高，計算解的精確性可能越好。的分析，可知的分析，可知Euler法具有一階精度，梯形法具二階精度。法具有一階精度，梯形法具二階精度。通過簡單通過簡單),(1nnnnutfhutu下面利用下面利用Taylor展開，求展開，求Euler法的局部截斷誤差法的局部截斷誤差

13、111)()(nnnutuhR nnnntutfhtutu,1 nnntuhtutu1 232!nnnhu thu tutO h nntuhtu 232!nhutO h2O h),(1nnnnuthfuu 線性二步法公式：線性二步法公式：)4(3122nnnnnfffhuu在前面，我們介紹了基于數(shù)值積分的特殊的單步法與多步法。在前面，我們介紹了基于數(shù)值積分的特殊的單步法與多步法。顯式單步法顯式單步法Euler公式公式: :單步法一般可寫成：單步法一般可寫成：111( ,; )nnnnnnuuht u tuh其中其中是依賴于是依賴于( (7-1) )右端的函數(shù)右端的函數(shù) f(t,u)。 當(dāng)取當(dāng)取

14、= =f( (tn,un) )時，為時，為Euler法法； 當(dāng)取當(dāng)取= =f( (tn+1,un+1) )時，為時，為隱式隱式Euler法法； 當(dāng)取當(dāng)取 時，為時，為梯形法梯形法。 111( ,)(,)2nnnnf tuf tu0,1, 2,1nN用到前一結(jié)點的值用到前一結(jié)點的值un-1-1, , 所以從初值所以從初值u0出發(fā)可逐步算出以后各結(jié)點出發(fā)可逐步算出以后各結(jié)點的值的值u1, ,u2, ,通過通過（7-3）計算結(jié)點計算結(jié)點tn=t0+nh, n=0,1,2,=0,1,2,的近似值的近似值 un , ,每次只每次只故稱為故稱為單步法單步法?；诨赥aylor展開的解法展開的解法線性多步

15、法的一般公式為：線性多步法的一般公式為：是常數(shù)，是常數(shù)， 0,00kkjjnkjjjnjfhu（7-4）其中其中 jjjnjnjnutff,),(0和和 0不同時為不同時為0 0。 按按（7-4）計算計算un+ +k時要用到前面時要用到前面k個結(jié)點值個結(jié)點值因此稱因此稱（7-4）為為多步法多步法 或或 k- -步法步法。 11,nnn ku uu 又因為又因為（7-4）關(guān)于關(guān)于 jnjnfu ,是線性的，所以稱為是線性的，所以稱為線性多步法線性多步法。為使多步法的計算能夠進行，除給定的初值為使多步法的計算能夠進行，除給定的初值u0 外，還要外，還要知道附加初值知道附加初值u1,u2, , ,u

16、k-1 ，這可用其它方法計算。，這可用其它方法計算。 若若k0,0,則方法（則方法（7-4）是）是隱式的隱式的。則稱（則稱（7-4）是）是顯式的顯式的；若若k=0例如，對于例如，對于線性二步法：線性二步法：nnnnnnfffhuuu0112201122當(dāng)取當(dāng)取 01,10,21,021,314,3時，就是時，就是Milne法。法。00kjnjjnjjuhf或或（數(shù)值解滿足的差分方程）（數(shù)值解滿足的差分方程）我們將基于我們將基于Taylor展開式來構(gòu)造出更一般的求解公式。展開式來構(gòu)造出更一般的求解公式。1、 基于基于Taylor展開式的求解公式展開式的求解公式 ( (單步法、多步法單步法、多步法

17、) )2、 顯式顯式Runge-Kutta法法( (單步法單步法) )用數(shù)值積分法只能構(gòu)造一類特殊的多步法，本節(jié)用數(shù)值積分法只能構(gòu)造一類特殊的多步法，本節(jié)這將包含如下的內(nèi)容：這將包含如下的內(nèi)容：初值問題初值問題（7-1），（），（7-2）的解充分光滑，將的解充分光滑，將 u(t) 在在t0處處用用Taylor公式展開：公式展開：10)(0200)(!)(! 2)()()( ppphOtuphtuhtuhtutu( (7-5) ) 其中其中,)(00utu),()(,()(00000utftutftu00( )t tdfu tdt0)(0) 3(ttdtdfdtdtu),(),(),(),(2)

18、,(00002000000utfutfutfutfutfuututt200000000),(),(),(),(utfutfutfutfuut( (7-6) ) 00( ,)tf t u00( ,)uf t u0( , )tt td f t udt0( , )( , )ut tdf t u f t udt000000( , )( , ) ( , )tttuf t uf t u f t u000000( ,)( ,)( ,)tuf t uf t uf t u00( ,)uf t u000000( ,)( ,) ( ,)utuuft uft uf t u00( ,)f t udtdt00( ,)f

19、t u dttdu0令令舍去余項舍去余項( ,( ); )t u th( (7-7) ) 則可將（則可將（7-5）改寫成為）改寫成為)();(,()()(10000phOhtuthtuhtu，則得，則得)(1phO1000( ,; )uuhtuh。一般而言，若已知一般而言，若已知un，則則1( ,; )nnnnuuhtuh 這是一個單步法，局部截斷誤差為這是一個單步法，局部截斷誤差為O( (h p+1+1) )， (7-7),(7-7),由（由（7-6） 當(dāng)當(dāng) p=1=1時，它是時，它是Euler法。法。展開法做數(shù)值計算，但可用它計算附加值。展開法做數(shù)值計算，但可用它計算附加值。 0,1, 2

20、,n 1pj( )jd u tjdt1jh! j11pj1( , ( )jdf t u t1jdt1jh! j1，是非線性的是非線性的關(guān)于關(guān)于可知可知f一般不直接用一般不直接用Taylor ，的工作量過大的工作量過大由于計算由于計算hutnn,令令( );kL u t h ( (7-28) ) 設(shè)設(shè)u( (t) )是初值問題是初值問題（7-1），（），（7-2）的解，將的解，將u(t+jh)和和 ()u tjh 在點在點t處進行處進行Taylor展開展開, ,()u tjh()u tjh 將上式代入（將上式代入（7-28）式，）式，得得( );kLu th 2.2.待定系數(shù)法待定系數(shù)法0kj(

21、)()jju tjhhu tjh23(3)()()( )( )( )( )1!2!3!jhjhjhu tu tu tut 23(3)(4)()()( )( )( )( )1!2!3!jhjhjhu tu tutut0kj( )0( )!llljhutl1( )1( )1 !llljhutljhj( )0( )!llljhutl1( )1( )1 !llljhutl0kj0kjjhj( )u t( )u t1!jh2()2!jh( )u t3()3!jh(3)( )ut( )u t( )ut1!jh22!j h( )ut33!j h(4)( )ut0kjj( )u t231!jh( )u t2(

22、)2!jh( )u t3()3!jh(3)( )ut234jjjh0kjjj! 220kjjj2h將下式按將下式按h的同次冪合并同類項，的同次冪合并同類項，得得0( )kjju t1l( )( )llut h記為記為( );kLu th 2()012 ( )( )( )( )pppc u tc hu tc h u tc h ut10011!1 !kklljjjjjjll10011!1 !kklljjjjjjll其中其中l(wèi)c, 1 , 0pl 11pphOc即有即有( (7-30) ) 若若u(t)有有p+2+2次連續(xù)微商，則可選取適當(dāng)次連續(xù)微商，則可選取適當(dāng)k和和j，j使使 0120 ,pcc

23、cc而而cp+1+10,0,即選即選j和和 j滿足：滿足：121(2)!pppkckp1122kck001kc2,3,p 01()k11121(2)1 !ppkkp121(2)!ppkkp122kk010k2,3,p 01()0k11121(2)01 !ppkkp此時此時)()( );(2) 1(11ppppkhOtuhchtuL而而 ，則，則)(,()(tutftunkjnnjnjRjhtujhtfhjhtu0)(,()()()( 2)1(11pnpppnhOtuhcR( );n kkRL u t hcp+1h p+1u(p+1)(t)稱為稱為局部截斷誤差主項局部截斷誤差主項; ; cp+1

24、稱為局部截斷誤差稱為局部截斷誤差主項系數(shù)主項系數(shù)。舍去余項舍去余項Rn，并，并un+j代替代替u(tn+jh)，用，用 fn+j記記f(tn+j,un+j)，就得到就得到 線性多步法線性多步法（7-24），），其其局部截斷誤差為局部截斷誤差為： 1(1)21 ( )()ppppchutO h可以證明其整體截斷誤差可以證明其整體截斷誤差En= =O( (h p) )，故稱為，故稱為p階階k步法步法。因為（因為（7-24）可以相差一個非零常數(shù)，所以不妨設(shè)）可以相差一個非零常數(shù)，所以不妨設(shè)k=1。 當(dāng)當(dāng)k=0時時, ,un+k可用可用un+k-1-1, ,，un直接表示，故稱為直接表示，故稱為顯式法

25、顯式法。 反之，當(dāng)反之，當(dāng)k00時，求時，求un+k需解一個方程需解一個方程( (一般用迭代法一般用迭代法),),用待定函數(shù)法構(gòu)造多步法的一個基本要求是選取用待定函數(shù)法構(gòu)造多步法的一個基本要求是選取j, j使使局部截斷誤差的階盡可能高。局部截斷誤差的階盡可能高。 下面我們討論構(gòu)造一般線性二步法公式的待定系數(shù)法。下面我們討論構(gòu)造一般線性二步法公式的待定系數(shù)法。22,1k 。記記0= =, ,其余四個系數(shù)其余四個系數(shù) 1012,由由 03210cccc確定，即滿足方程：確定，即滿足方程：稱為稱為隱式法隱式法。此時此時nnnnnnfffhuuu0112201122解之得解之得 所以一般二步法為所以一

26、般二步法為21(1)nnnuuu( (7-31) ) 0110c 110122()0c21121(4)(2)02c311211(8)(4)062c1(1), 01(1 5 ),12 21(5),1212(1),321(5)8(1)(1 5 )12nnnhfff121(2)!pppkckp11121(2)1 !ppkkp由（由（7-30）知道）知道: :411211(16)(8)246c511211(32)(16)12024c當(dāng)當(dāng)-1-1時時c40 ，方法（，方法（7-31）是三階二步法。）是三階二步法。 )4(3122nnnnnfffhuu這是四階二步法，是具有最高階的二步法，稱為這是四階二步

27、法，是具有最高階的二步法，稱為Milne方法方法。nnnnnfffhuu121285121(1) ,24 1(17 13 ),360 當(dāng)當(dāng)=-1=-1時時c4= =0 ，但，但c50，方法（方法（7-31）化為：）化為： 此外，若取此外，若取= =0, ,則（則（7-31）為：）為： 此為二步隱式此為二步隱式Adams方法；若取方法；若取=-5 =-5 ，則（，則（7-31）為：）為： 2114522nnnnnuuuhff是顯式方法。是顯式方法。用類似的計算過程可獲得一些常用的線性多步法的局部用類似的計算過程可獲得一些常用的線性多步法的局部截斷誤差。截斷誤差。當(dāng)當(dāng)k=1=1時，梯形法（二階隱式方法）：時，梯形法（二階隱式方法）： nnnnffhuu112其中其中1=1, , 0=-