chap常微分方程數(shù)值解PPT課件

1、第七章第七章 常微分方程數(shù)值解法常微分方程數(shù)值解法1 、引言第1頁/共120頁2yxy2y4xy4x2yf x y4y 13xdy2yf x y4dxxy 13: - ( , ) ( )-( , )( ) 例例方方程程令令：且且給給出出初初值值就就得得到到一一階階常常微微分分方方程程的的初初值值問問題題：第2頁/共120頁 ( , )( , )( , )f x yyLipschitzLf x yf x yL yy只要函數(shù)適當光滑連續(xù)，且關于 滿足條件，即存在常數(shù) ，使得由常微分方程理論知，初值問題的解必存在且唯一。第3頁/共120頁 微分方程的數(shù)值解：設方程問題的解y(x)的存在區(qū)間是a,b，

2、令a= x0 x1 xn =b,其中hk=xk+1-xk , 如是等距節(jié)點h=(b-a)/n , h稱為步長。y(x)的解析表達式不容易得到或根本無法得到，我們用數(shù)值方法求得y(x)在每個節(jié)點xk上y(xk)的近似值，用yk表示，即yky(xk)，這樣y0 , y1 ,.,yn稱為微分方程的數(shù)值解。 第4頁/共120頁主要問題v如何將微分方程離散化，并建立求其數(shù)值解的遞推公式；v遞推公式的局部截斷誤差，數(shù)值解與精確解的誤差估計；v遞推公式的穩(wěn)定性與收斂性。第5頁/共120頁v用差商代替微商用差商代替微商v數(shù)值積分數(shù)值積分vTaylor展開展開微分方程離散化常用方法微分方程離散化常用方法第6頁/

3、共120頁第7頁/共120頁第8頁/共120頁第9頁/共120頁第10頁/共120頁第11頁/共120頁1 解常微分方程初值問題的Euler方法Euler方法Euler方法的誤差分析第12頁/共120頁v向前向前EulerEuler公式公式（EulerEuler折線法或顯格式）折線法或顯格式）v向后Euler公式（后退Euler公式）v梯形公式（改進的Euler公式）vEuler預估校正格式一、Euler方法第13頁/共120頁1、向前Euler公式第14頁/共120頁第15頁/共120頁用分段的折線逼近逼近函數(shù)第16頁/共120頁2、向后（后退的）Euler 方法第17頁/共120頁第18頁

4、/共120頁3、梯形公式第19頁/共120頁第20頁/共120頁4、改進的尤拉公式梯形公式雖然提高了精度，但使算法復雜。而在實際計算中只迭代一次，這樣建立的預測校正系統(tǒng)稱作改進的尤拉公式。第21頁/共120頁1111 ( ,); ( ,)(,),2nnnnnnnnnnyyhf x yhyyf x yf xy預測校正11(,);(,);() / 2.pnnncnnpnpcyyhfxyyyhfxyyyy第22頁/共120頁二、二、Euler方法的誤差分析方法的誤差分析第23頁/共120頁第24頁/共120頁第25頁/共120頁第26頁/共120頁2） 總體方法誤差總體方法誤差第27頁/共120頁第

5、28頁/共120頁第29頁/共120頁第30頁/共120頁總體截斷誤差與局部截斷誤差的關系是：第31頁/共120頁誤差分析表誤差分析表EulerEuler方法方法局部截局部截斷誤差斷誤差總體截總體截斷誤差斷誤差迭代收斂迭代收斂條件條件向前向前EulerEuler方法方法O(h2)O(h)向后向后EulerEuler方法方法O(h2)O(h)0hL1梯形公式梯形公式O(h3)O(h2)0hL2(L為為Lip常數(shù)）常數(shù)）第32頁/共120頁向后Euler 方法收斂條件與截斷誤差第33頁/共120頁梯形公式的收斂性第34頁/共120頁2(01);(0)1.xyyxyy例：用尤拉公式和改進的尤拉公式解

6、初值問題10.12().nnnnnhxEuleryyh yy解：取步長，公式為： 第35頁/共120頁112();2();1().2npnnnncnppnpcxyyh yyxyyh yyyyy改進的尤拉公式為：第36頁/共120頁20.2(00.6); (0)1.hyyxyxy 例1：取步長，用歐拉法解初值問題122 , 0.2 0.80.2nnnnnnnnnnnEuleryyhfxyyyx yyx y解：格式為：第37頁/共120頁0122 1 0.2 0.8, 0.4 0.6144 0.6 0.461321yyyyyyy由計 算 得第38頁/共120頁0.283(12); (1)2.5hy

7、yxy例2：取步長，用梯形解初值問題小數(shù)點后至少保留 位。11111 ,20.2 83832nnnnnnnnnnhyyfxyfxyyyyy解：梯形公式為： 第39頁/共120頁 1012345716 131312, 1.22.307691.42.473371.62.562581.82.610622.02.63649nnyyyyyyyyyyyyyy故 由計 算 得 第40頁/共120頁0.23(12); (1)2.5hyxyxy例2 ：取步長，用梯形解初值問題小數(shù)點后至少保留 位。111111 ,2 0.1 33nnnnnnnnnnnnhyyfxyfxyyyx yxy解：梯形公式為：第41頁/共

8、120頁 0111110011234511111110.30.1 332,2.6*,&,%,?nnnnkknnnnnnkyyx yyyx yxyyyyyyyyyy迭代格式：第42頁/共120頁 021234522222,*,&,%,yyyyyy第43頁/共120頁20.2sin0 (1)1.1.21.45hyyyxyyy例3：取步長，用歐拉預校方法解初值問題計算及的近似值，小數(shù)點后至少保留 位。第44頁/共120頁111121212111,2 0.2sin0.1sin + sinnnnnnnnnnnnnnnnnnnnnnnnyyhfxyhyyfxyfxyyyyyxyyyyxyy

9、x解：歐拉預校格式為：第45頁/共120頁 0112211,0.63171 1.20.7154880.476961.40.52611yyyyyyyy故 由計 算 得 第46頁/共120頁50 (0)1.0.10.210yyyhy例4：寫出用反復迭代的歐拉預校法解初值問題的計算公式，并取步長，計算，要求迭代誤差不超過。第47頁/共120頁 01111101111,2 0,1,2, 0,1,2, 0.90.950.05nnnnkknnnnnnnnkknnnyyhfxyhyyfxyfxyknfx yyyyyyy 解：歐拉預校的迭代格式為：取第48頁/共120頁 0012111341143751141

10、101, 0.9,0.905,0.90475, 0.9047625,0.9047618756.251010, 0.10.904761875yyyyyyyyyyyy故 由計 算 得由于 是 取第49頁/共120頁 0012232452254652252201, 0.814286,0.818809, 0.818583, 0.818595,0.8185941010, 0.20.818594yyyyyyyyyyyy故 由計 算 得由于 是 取第50頁/共120頁2.2.龍格龍格庫塔方法庫塔方法基本思想基本思想二階二階R-K方法方法三階三階R-K方法方法四階四階R-K方法方法變步長變步長R-K方法方法第

11、51頁/共120頁1112() ()()()()()()2!( )( , ), ( )( , )( , ) ( , ),nnnppnnnnxypTaylory xyy xhhy xhy xyxyxPy xf x yyxfx yfx y f x y若用 階多項式近似函數(shù)有：其中。但由于公式中各階偏導數(shù)計算復雜，不實用。一、基本思想第52頁/共120頁(0)(0)(0)(1)(1)(1)(2)(2)(2)( )(1); ;2,3,jjjjyffffyffxyffyffxyffyffjxy一般地有第53頁/共120頁111112121(,)11()22 (,)(,)nnnnnnnnnnEulerEu

12、lerEuleryyhKKf xyEuleryyhKKKf xyKf xh yhK如果將公式與改進公式寫成下列形式：公式改進公式第54頁/共120頁11 ( , )()( , )( , )nnf x yy xyf x yf x y以上兩組公式都使用函數(shù)在某些點上的值的線性組合來計算的近似值。Euler公式：每步計算一次的值，為一階方法。改進Euler公式：需計算兩次的值，二階方法。第55頁/共120頁( , )(,)( )-nnnf x yxyTaylory xxTaylorR K于是可考慮用函數(shù)在若干點上的函數(shù)值的線性組合來構造近似公式，構造是要求近似公式在處的展開式與解在 處的展開式的前面

13、幾項重合，從而使近似公式達到所需要的階數(shù)。即避免求偏導，又提高了方法的精度，此為方法的基本思想。第56頁/共120頁11111-(,)(,) (2,3,),pnniiinniininijjjiijiRKyyhc KKf xyKf xa h yhb Kipabc一般地，方法設近似公式為其中 ，都是參數(shù).第57頁/共120頁,(,)( )iijinnnab cxyTaylory xxTaylor確定參數(shù) ，的原則是:使近似公式在處的展開式與在 處的展開式的前面項盡可能多地重合。第58頁/共120頁二、二階龍格庫塔方法111221222112()(,)(,)nnnnnnpyyh c Kc KKf x

14、yKf xa h yhb K當時，近似公式為 第59頁/共120頁112221123221( ,)( ,)(, ( ,)( ,) ( ,) ( ,)( ,) ( ,)( )nnnnnnnnnnnnnnnxnnynnnnnx yTayloryyhc f x yc f xa h yhb f x yyh c f x yc f x ya hf x yhb f x y f x yO hy上式在處的展開式為12222321() ( ,)( ,) + ( ,) ( ,)( )nnxnnynnnncc f x y h c a f x yb f x y f x yhO h第60頁/共120頁123123()()

15、()()()()2(,)(,)2 (,) (,)()nnnnnnnnnxnnynnnny xxTaylorhy xy xhy xy xO hhyf xy hfxyfxyf xyO h在 處的展開式為第61頁/共120頁122 232 211 1/2 1/2 ( ),cccacbOh有 無 窮 多 組 解 ， 每 一 組 解 得 一近 似 公 式 ， 局 部 截 斷 誤 差 均 為這 些 方 法 統(tǒng) 稱 二 階 方 法 。43()()O hO h可以證明，無論這四個參數(shù)如何選擇，都不能使局部截斷誤差達到，也即在計算兩次函數(shù)值的情況下，局部截斷誤差的階最高為。第62頁/共120頁122211121

16、211,1,2()/2(,)(,)nnnnnnccabEuleryyh KKKf xyKf xh yhK取此為改進公式。近似公式為 第63頁/共120頁122211212110,1,2( ,)(2,2)nnnnnnccabyyhKKf x yKf xhyhK取此為常用的二階公式，稱為中點公式 第64頁/共120頁三、三階龍格庫塔方法1123121312(4)6(,)(,)22(,2)nnnnnnnnRKhyyKKKKf xyhhKf xyKKf xh yhKhK常用的三階公式為：第65頁/共120頁四、四階龍格庫塔方法112341213243 (22)6(,)(,)22(,)22(,)nnnn

17、nnnnnnRKhyyKKKKKf xyhhKf xyKhhKf xyKKf xh yhK常用的四階公式為： 第66頁/共120頁0.2,-83 ;(0)2.0.44 hR Kyyyy例：設取步長寫出用經(jīng)典（標準的）四階方法求解初值問題 的計算公式，計算的近似值，小數(shù)點后至少保留 位。12,8 - 3 , =0.2, 0.2,0.4fx yyhyyyy解：第67頁/共120頁112341122343(22);6,83;,5.62.1;22,6.322.37;22,4.2081.578.nnnnnnnnnnnnnnhyyKKKKKf xyyKhKfxyyKhKfxyyKf xh yKy由經(jīng)典的四

18、階龍格庫塔公式得第68頁/共120頁 10121.20160.549402,0.22.30040.42.4654nnyyyyyyyy由于第69頁/共120頁11,2,3,4,454652pR KppR KpppR KR KTaylor兩點說明：）當時， 公式的最高階數(shù)恰好是 當時， 公式的最高階數(shù)不是，如時仍為 ，時 公式的最高階數(shù)為。） 方法的導出基于展開，故要求所求問題的解具有較高的光滑度。第70頁/共120頁 RKEulerRKRKEuler當解充分光滑時，四階 方法確實優(yōu)于改進法。對一般實際問題，四階方法一般可達到精度要求。如果解的光滑性差，則用四階 方法解的效果不如改進法。第71頁/

19、共120頁五、變步長的龍格庫塔方法( )1( )51115(2)1,(),2,2nhnhnnnnhnxhyy xychhxxhyc 以經(jīng)典四階龍格庫塔公式為例。從節(jié)點 出發(fā)，以 為步長求一近似值將步長折半，即取 為步長從 跨兩步到，求一近似值每跨一步的截斷誤差是第72頁/共120頁5(2)11(2)11()11(2)(2)()1111()2,2()1.16()1().15hnnhnnhnnhhhnnnnhy xycy xyy xyy xyyy因此有由上兩式 第73頁/共120頁 4、微分方程數(shù)值解的穩(wěn)定性第74頁/共120頁第75頁/共120頁Euler法的絕對穩(wěn)定區(qū)域第76頁/共120頁 I

20、m(h)-2 -1 0 Re(h)第77頁/共120頁向后向后Euler Euler 法的穩(wěn)定性法的穩(wěn)定性第78頁/共120頁第79頁/共120頁梯形公式的穩(wěn)定性第80頁/共120頁第81頁/共120頁R-K方法的絕對穩(wěn)定區(qū)域第82頁/共120頁第83頁/共120頁 2 1-3 -2 -1 0 -1 -2第84頁/共120頁11-r1R-K,nnnnnnyyyyyy單步法在計算時，只用到前一步的信息 。為提高精度，需重新計算多個點處的函數(shù)值，如方法，計算量較大。如何通過較多地利用前面的已知信息，如 ，來構造高精度的算法計算?；舅枷?.線性多步法第85頁/共120頁11110111 (,),

21、(,) ,( ,)(1, ,)00Taylornnn rnnn rn rrrnin iin iiiiikkkyyyf xyf xyyyhfff x yknnn r 多步法中最常用的是線性多步法，它的計算公式中只出現(xiàn)， ，及的一次項，其一般形式為 其中均為常數(shù)，。若，顯式；，隱式。構造線性多步公式常用展開和數(shù)值積分方法。第86頁/共120頁線性多步公式的導出1(),nnniiTaylorxTaylory xxTaylor 利用展開導出的基本方法是：將線性多步公式在 處進行展開，然后與在處的展開式相比較，要求它們前面的項重合，由此確定參數(shù)。第87頁/共120頁( )( )2( )1() (1,2,

22、),( )( )()()2 ()()!kknnnnnnnnpppnnnyyxky xxTayloryy xyy xxxxyxxO xxp記則在 處的展開為1011110111( ) ()nnnnnnry xyyyhfff以為例：設初值問題的解充分光滑，待定的兩步公式為第88頁/共120頁231(4)(5)45(6)21111(),()( ,) (),()2!3! ()4!5!(,)()2! iiiiinnnnnnnnnnnnnnnnyy xy xf x yinyyyy xhyy hhhyyhhO hyff xyy xyy hh假設前 步計算結果都是準確的，即則有(4)(5)34(5) ()3!

23、4!nnyyhhO h第89頁/共120頁1111(4)(5)234(5)(,)(,) () ()2!3!4!nnnnnnnnnnnnnff xyyff xyy xyyyyy hhhhO h第90頁/共120頁1(5)2561()()() 2!5!1 nnnnnnnpy xxTayloryyy xyy hhhO hp為使上式有 階精度，只須使其與在 處的展開式的前項重合。211011101113(4)4111111(5)56111()()()2()()6222466()( )1202424nnnnnnnyyy hy hy hy hy hO h 將以上各公式代入并整理，得第91頁/共120頁01

24、0101111111111111122111162261111246624aaaaaa5,5,11,2,3,4)iippp 個參數(shù)只須 個條件。由推導知，如果選取參數(shù)，使其滿足前個方程（，則近似公式為 階公式。第92頁/共120頁11011111,0,02 ()2nnnnhyyff0如滿足方程組前三個方程，故公式此為二階公式。第93頁/共120頁0111011115(5)6110,1,343 (4)31 ( )90nnnnnnnhyyfffRh yOh又如：解上面方程組得,相應的線性二步四階公式（Simpson公式)為其截斷誤差為由此可知，線性二步公式至多是四階公式。第94頁/共120頁101

25、0123()1( )( )1 (1,2,)nnnriirrkkiiiirxTaylory xxTaylorikikp一般地，線性多步公式中有個待定參數(shù)，如令其右端在 處的展開式與在處的展開式的前p+1項系數(shù)對應相等，可得方程組第95頁/共120頁1111(1)21 1( )(1)! (1)( ) ()prpniirpppiniphRippiyO h其解所對應的公式具有 階精度，局部截斷誤差為顯然，線性多步公式至多可達到2r+2階精度。第96頁/共120頁二、常用的線性多步公式第97頁/共120頁1231010100123 (Adams)r=30,1()()1 (1,2,3 4)5559379=

26、1,24242424riirrkkiiiiikik （一）阿達姆斯公式取，并令由方程組，可解得，第98頁/共120頁1123153354(5)61115(5)6(5559379)24=0Adams1( )5( )()5!251()720nnnnnnniiniinhyyffffhRiiyO hh yO h相應的線性多步公式為因，此式稱為顯式公式，是四階公式.局部截斷誤差為第99頁/共120頁12330101211125(5)610,91951 =1,24242424(9195)24Adams19()720nnnnnnnnhyyffffRh yO h如果令由方程組可解得，相應的線性多步公式為稱其為

27、四階隱式公式，其局部截斷誤差為第100頁/共120頁利用數(shù)值積分方法求線性多步公式111111()()( , ( )( ),( )( )( ),1nnnnxxnnxxnnn rnnn rry xy xf x y x dxF x dxxxxxxxF xrxF xr 基本思想是首先將初值問題化成等價的積分形式用過節(jié)點或的的 次插值多項式代替求積分即得階的線性多步公式。第101頁/共120頁123330123303,( ) ( )( ) ()()()()()( )()() (0,1,2,3)nnnnin iinnnnin in injjj irx xxxF xL xl x F xxxxxxxxxl

28、xxxxxi例如時，過節(jié)點的三次插值多項式為其中第102頁/共120頁1111131301233231313233()()( )( ) ()()()()()6()()()()2()()()()2()()nnnnnnnnnnxxnnin ixxixnnnnxxnnnnxxnnnnxnny xy xL x dxl x dx F xxxxxxxF xdxhxxxxxxF xdxhxxxxxxF xdxhxxxF x1123123)()655 ()59 ()37 ()9 ()24nnxnnxnnnnxxxdxhhF xF xF xF x第103頁/共120頁1111233,(), (),(,)()(,

29、 () (,1,2,3),(5559379)24,nnnnkkkkkknnnnnnnnyyy xy xfxyF xf xy xkn nnnhyyffffAdamsxxAdams對上式用代替用代替則得這就是四階顯式公式。由于積分區(qū)間在插值區(qū)間外面，又稱為四階外插公式。第104頁/共120頁111(4)310(5)3031(5)35(5)10()()4!() ()4!,),( )251()( )4!720nnnnnnxxnnjxjxxnjxjnnxnnjxjFRxxdxyxxdxxxyRxxdxh y由插值余項公式可得其局部截斷誤差為由積分中值定理，存在（使得第105頁/共120頁11223111

30、231,( ) ( )( ) ()()()()()( )()() (1,0,1,2)nnnnin iinnnnin in in jjj ixx xxF xL xl x F xxxxxxxxxl xxxxxi 同樣，如果過節(jié)點的三次插值多項式為其中代( )F x替求積分，第106頁/共120頁11125(5)12121 (9195)2419( )720,nnnnnnnnnnnnAdamshyyffffRh yxxxxAdamsAdams 即得四階隱式公式其局部截斷誤差為由于積分區(qū)間在插值區(qū)間內(nèi)，故隱式公式又稱為內(nèi)插公式第107頁/共120頁一階常微分方程組與高階方程 我們已介紹了一階常微分方程初

31、值問題的各種數(shù)值解法，對于一階常微分方程組，可類似得到各種解法，而高階常微分方程可轉化為一階常微分方程組來求解。一階常微分方程組對于一階常微分方程組的初值問題 0000( , , ), ()( , , ), ()yf x y z y xyzg x y z z xz（） 可以把單個方程 中的f 和y看作向量來處理，這樣就可把前面介紹的各種差分算法推廣到求一階方程組初值問題中來。 ( , )yf x y 第108頁/共120頁設 為節(jié)點上的近似解，則有改進的Euler格式為 0(1,2,3,);,iiixxih iy z1( ,)iiiiiyyhf x y z1( ,)iiiiizzhg x y

32、z預報：1111( ,)(,)2iiiiiiiihyyf x y zf xyz1111( ,)(,)2iiiiiiiihzzg x y zg xyz校正： （） 又，相應的四階龍格庫塔格式（經(jīng)典格式）為 1123411234(22)6(22)6iiiihyyKKKKhzzLLLL（） 第109頁/共120頁112111221112312223122241334133(,)(,)(,)22(,)22(,)22(,)22(,)(,)iiiiiiiiiiiiiiiiiiiiiiiiKfxyzLg xyzhhKfxyKzLhhLg xyKzLhhKfxyKzLhhLg xyKzLKfxyhKzhLLg

33、 xyhKzhL式中 （） 第110頁/共120頁 把節(jié)點xi上的yi和zi值代入式(7.34), 依次算出 , 然后把它們代入式(7.33), 算出節(jié)點xi+1上的yi+1 和zi+1值。 對于具有三個或三個以上方程的方程組的初值問題,也可用類似方法處理,只是更復雜一些而已。此外,多步方法也同樣可以應用于求解方程組初值問題。 11223344,K L KL KL KL例 用改進的Euler法求解初值問題 (0)1(0)2yxyzyxyzzz 00.2x取步長，保留六位小數(shù)。 解: 改進的Euler法公式為),(1iiiiizyxhfyy),(1iiiiizyxhgzz預報： ),(),(21

34、111iiiiiiiizyxfzyxfhyy),(),(21111iiiiiiiizyxgzyxghzz校正： 將 及代入上式,得 zyxzyxgzxyzyxfiiiiii),(,),(第111頁/共120頁110.1()0.1iiiiiiiiiiyyx yzxyzzz11111110.05 ()()0.05iiiiiiiiiiiiiiiiyyx yzxyzxyxyzzzz由初值 ,計算得 00(0)1,(0)2yyzz110.8000002.050000yz11(0.1)0.801500(0.1)2.046951yyzz220.6048202.090992yz22(0.2)0.604659(0.2)2.088216yyzz第112頁/共120頁高階方程組 高階微分方程(或方程組)的初值問題,原則上都可以歸結為一階方程組來求解。例如,有二階微分方程的初值問題 0000( , ,)(),()yf x y yy xyy xy在引入新的變量 后,即化為一階方程組初值問題:zy0000( , , ), (), ()zf x y zyz y xy z xy 式（）為一個一階方程組的初值問題，對此可用中介紹的方法來求解。例如應用四階龍格-庫塔