1、工程數(shù)學(xué)工程數(shù)學(xué)工程數(shù)學(xué)工程數(shù)學(xué)1第十章第十章 常微分方程數(shù)值解常微分方程數(shù)值解第一節(jié)第一節(jié) 求解初值問題數(shù)值方法的基本原理求解初值問題數(shù)值方法的基本原理第二節(jié)第二節(jié) 高精度的單步法高精度的單步法 第三節(jié)第三節(jié) 線性多步法線性多步法第四節(jié)第四節(jié) 一階微分方程組的解法一階微分方程組的解法第五節(jié)第五節(jié) 邊值問題的打靶法和差分法邊值問題的打靶法和差分法工程數(shù)學(xué)工程數(shù)學(xué)工程數(shù)學(xué)工程數(shù)學(xué)2考慮一階常微分方程的初值問題考慮一階常微分方程的初值問題 /* Initial-Value Problem */: 0)(,),(yaybaxyxfdxdy只要只要 f (x, y) 在在a, b R1 上連續(xù)，且關(guān)于

2、上連續(xù)，且關(guān)于 y 滿足滿足 Lipschitz 條條件件，即存在與，即存在與 x, y 無關(guān)的常數(shù)無關(guān)的常數(shù) L 使使對任意定義在對任意定義在 a, b 上的上的 y1(x) 和和 y2(x) 都成立，則上述都成立，則上述IVP存存在唯一解。在唯一解。| ),(),(|2121yyLyxfyxf 要計算出解函數(shù)要計算出解函數(shù) y(x) 在一系列節(jié)點在一系列節(jié)點 a = x0 x10,使得使得),(hyxhQyynnnn 111()()pO hp (, )nnQ xyh( , , )( , , )Q x y hQ x y hL yy對一切對一切 成立成立,則該方法收斂則該方法收斂,且有且有 y

3、y和和)(pnhOe 工程數(shù)學(xué)工程數(shù)學(xué)工程數(shù)學(xué)工程數(shù)學(xué)14由該定理可知整體截斷誤差總比局部截斷誤差低一階由該定理可知整體截斷誤差總比局部截斷誤差低一階 對改進的對改進的Euler法法, ),(,(),(),(yxhfyhxfyxfhyxQ 21于是有于是有 1( , , )( , , )( , )( , )2(,( , )(,( , )Q x y hQ x y hf x yf x yf xh yhf x yf xh yhf x y設(shè)設(shè)L為為f關(guān)于關(guān)于y的的Lipschitz常數(shù)常數(shù),則由上式可得則由上式可得( , , )( , , )(1/2)Q x y hQ x y hLhyy限定限定h即可

4、知即可知Q滿足滿足Lipschitz條件條件,故而改進的故而改進的Euler法收斂法收斂.工程數(shù)學(xué)工程數(shù)學(xué)工程數(shù)學(xué)工程數(shù)學(xué)15例：例：考察初值問題考察初值問題 在區(qū)間在區(qū)間0, 0.5上的解。上的解。分別用歐拉顯、隱式格式和改進的歐拉格式計算數(shù)值解。分別用歐拉顯、隱式格式和改進的歐拉格式計算數(shù)值解。 1)0()(30)(yxyxy0.00.10.20.30.40.5精確解精確解改進歐拉法改進歐拉法 歐拉隱式歐拉隱式歐拉顯式歐拉顯式 節(jié)點節(jié)點 xixey30 1.0000 2.0000 4.0000 8.0000 1.6000 101 3.2000 101 1.00002.5000 10 1 6

5、.2500 10 21.5625 10 23.9063 10 39.7656 10 41.00002.50006.25001.5626 1013.9063 1019.7656 1011.00004.9787 10 22.4788 10 31.2341 10 46.1442 10 63.0590 10 73. 穩(wěn)定性穩(wěn)定性工程數(shù)學(xué)工程數(shù)學(xué)工程數(shù)學(xué)工程數(shù)學(xué)16定義定義若某算法在計算過程中任一步產(chǎn)生的誤差在以后的計若某算法在計算過程中任一步產(chǎn)生的誤差在以后的計算中都算中都逐步衰減逐步衰減，則稱該算法是，則稱該算法是絕對穩(wěn)定的絕對穩(wěn)定的 /*absolutely stable */。一般分析時為簡單起

6、見，只考慮一般分析時為簡單起見，只考慮試驗方程試驗方程 /* test equation */yy 常數(shù)，可以常數(shù)，可以是復(fù)數(shù)是復(fù)數(shù)當步長取為當步長取為 h 時，將某算法應(yīng)用于上式，并假設(shè)只在初值時，將某算法應(yīng)用于上式，并假設(shè)只在初值產(chǎn)生誤差產(chǎn)生誤差 ，則若此誤差以后逐步衰減，就稱該，則若此誤差以后逐步衰減，就稱該算法相對于算法相對于 絕對穩(wěn)定絕對穩(wěn)定， 的全體構(gòu)成的全體構(gòu)成絕對穩(wěn)定區(qū)域絕對穩(wěn)定區(qū)域。我們稱我們稱算法算法A 比算法比算法B 穩(wěn)定穩(wěn)定，就是指，就是指 A 的絕對穩(wěn)定區(qū)域比的絕對穩(wěn)定區(qū)域比 B 的的大大。000yy h h h工程數(shù)學(xué)工程數(shù)學(xué)工程數(shù)學(xué)工程數(shù)學(xué)17例：例：考察顯式歐拉

7、法考察顯式歐拉法110(1)nnnnyyhyhy 000yy 110(1)nnyhy 11110(1)nnnnyyh 由此可見，要保證初始誤差由此可見，要保證初始誤差 0 以后逐步衰減，以后逐步衰減，必須滿足：必須滿足：hh 1|1| h0-1-2ReImg例：例：考察隱式歐拉法考察隱式歐拉法11nnnyyh y 111nnyyh 11011nnh 可見絕對穩(wěn)定區(qū)域為：可見絕對穩(wěn)定區(qū)域為：1|1| h210ReImg注：注：一般來說，隱式歐拉法的絕對穩(wěn)定性比同階的顯式一般來說，隱式歐拉法的絕對穩(wěn)定性比同階的顯式法的好。法的好。工程數(shù)學(xué)工程數(shù)學(xué)工程數(shù)學(xué)工程數(shù)學(xué)1812()11()()()()()

8、( )2!( )( ,),( )( ,)( ,)(), ,nppnnnnnnxypTaylory xhhyy xy xhyxyxyxPyxf x yyxfx yfx y f x y 若若用用 階階多多項項式式近近似似函函數(shù)數(shù)有有：其其中中。但但由由于于公公式式中中各各階階偏偏導(dǎo)導(dǎo)數(shù)數(shù)計計算算復(fù)復(fù)雜雜，不不實實用用。第二節(jié)第二節(jié) 高精度的單步法高精度的單步法在高精度的單步法中在高精度的單步法中, ,應(yīng)用最廣泛的是應(yīng)用最廣泛的是Runge-KuttaRunge-Kutta( (龍格龍格- -庫塔庫塔) )方法方法一一、Runge-Kutta法的基本思想（法的基本思想（1）工程數(shù)學(xué)工程數(shù)學(xué)工程數(shù)學(xué)工

9、程數(shù)學(xué)19(0)(0)(0)(1)(1)(1)(2)(2)(2)( )(1); ;2,3,jjjjyffffyffxyffyffxyffyffjxy 一般地有一般地有工程數(shù)學(xué)工程數(shù)學(xué)工程數(shù)學(xué)工程數(shù)學(xué)20111112121 (,)11()22(,)(, )nnnnnnnnnnE ulerE uleryyhKE ulerKfxyyyhKKE ulerKfxyKfxhyhK 如如 果果 將將公公 式式 與與 改改 進進公公 式式 寫寫 成成 下下 列列 形形 式式 ：公公 式式 改改 進進公公 式式Runge-Kutta法的基本思想（法的基本思想（2）工程數(shù)學(xué)工程數(shù)學(xué)工程數(shù)學(xué)工程數(shù)學(xué)2111( ,

10、)()( , )( , ) nnf x yy xyf x yf x y以上兩組公式都使用函數(shù)在某些點上的以上兩組公式都使用函數(shù)在某些點上的值的線性組合來計算的近似值。值的線性組合來計算的近似值。EulerEuler公式：每步計算一次的值，為一階方法。公式：每步計算一次的值，為一階方法。改進Euler公式：需計算兩次的值，二階方法。改進Euler公式：需計算兩次的值，二階方法。工程數(shù)學(xué)工程數(shù)學(xué)工程數(shù)學(xué)工程數(shù)學(xué)22Runge-Kutta法的基本思想（法的基本思想（3）( ,)(,)( ) nnnf x yxyTaylory xxTaylor于于是是可可考考慮慮用用函函數(shù)數(shù)在在若若干干點點上上的的函

11、函數(shù)數(shù)值值的的線線性性組組合合來來構(gòu)構(gòu)造造近近似似公公式式，構(gòu)構(gòu)造造是是要要求求近近似似公公式式在在處處的的展展開開式式與與解解在在處處的的展展開開式式的的前前面面幾幾項項重重合合，從從而而使使近近似似公公式式達達到到所所需需要要的的階階數(shù)數(shù)。即即避避免免求求偏偏導(dǎo)導(dǎo)，又又提提高高了了方方法法的的精精度度，此此為為R RK K方方法法的的基基本本思思想想。11111(,)(,)(2, 3,)pnniiinniininijjjyyhc KKfxyKfxa h yhb Kip 工程數(shù)學(xué)工程數(shù)學(xué)工程數(shù)學(xué)工程數(shù)學(xué)23二、二階龍格庫塔方法二、二階龍格庫塔方法11111RK(,)(,)(2,3,),(,)

12、()pnniiinniininijjjiijinnnyyhc KKfxyKfxa h yhb KipabcxyTaylory xxTaylor 一一般般地地，方方法法設(shè)設(shè)近近似似公公式式為為其其中中，都都是是參參數(shù)數(shù)，確確定定它它們們的的原原則則是是使使近近似似公公式式在在處處的的展展開開式式與與在在處處的的展展開開式式的的前前面面項項盡盡可可能能多多地地重重合合。工程數(shù)學(xué)工程數(shù)學(xué)工程數(shù)學(xué)工程數(shù)學(xué)2411122122211()2(,)(,)nnnnnnyyh c Kc KKf xyKf xa h yhb K 當p時，近似公式為 當p時，近似公式為 112221122321(,)(,)(,(,)

13、(,) (,)(,)(,) (,)()nnnnnnnnnnnnnnnxnnynnnnxyTayloryyh c f xyc f xa h yhb f xyyh c f xycf xya hfxyhb fxyf xyO h 上式在處的展開式為上式在處的展開式為工程數(shù)學(xué)工程數(shù)學(xué)工程數(shù)學(xué)工程數(shù)學(xué)2512232221()(,)(,)(,)(,)()nnnxnnynnnnyccf xyhc a fxybfxyf xyhO h 123123()()()()()()2(,)(,)(,)(,)()2nnnnnnnnnxnnynnnny xxTaylorhy xy xhyxyxO hyfxyhhfxyfxyfx

14、yO h 在在處處的的展展開開式式為為工程數(shù)學(xué)工程數(shù)學(xué)工程數(shù)學(xué)工程數(shù)學(xué)2612222211 1 / 2 1 / 2 ),ccc ac bO 3 3有有 無無 窮窮 多多 組組 解解 ， 每每 一一 組組 解解 得得 一一近近 似似 公公 式式 ， 局局 部部 截截 斷斷 誤誤 差差 均均 為為( (h h這這 些些 方方 法法 統(tǒng)統(tǒng) 稱稱 二二 階階 方方 法法 。122211121211,1,2() / 2(,)(,)nnnnnnccabE uleryyh KKKfxyKfxhyhK 取取此此 為為 改改 進進公公 式式 。近近 似似 公公 式式 為為 122211212110,1,2(,)

15、(2,2)nnnnnnccabyyhKKfxyKfxhyhK 取取此此 為為 常常 用用 的的 二二 階階 公公 式式 ，稱稱 為為 中中 點點 公公 式式 。 工程數(shù)學(xué)工程數(shù)學(xué)工程數(shù)學(xué)工程數(shù)學(xué)27三、三階龍格庫塔方法三、三階龍格庫塔方法11231213123 (4)6(,)( ,)22(,2)nnnnnnnnpRKRKhyyKKKKf xyhhKf xyKKf xh yhKhK 類類似似地地，對對，即即三三個個點點，通通過過更更復(fù)復(fù)雜雜的的計計算算，可可導(dǎo)導(dǎo)出出三三階階公公式式。常常用用的的三三階階公公式式為為：工程數(shù)學(xué)工程數(shù)學(xué)工程數(shù)學(xué)工程數(shù)學(xué)28四、四階龍格庫塔方法四、四階龍格庫塔方法11

16、23412132434 (22)6(,)(,)2 2(,)22(,)nnnnnnnnnnpRKRKhyyKKKKKf xyhhKf xyKhhKf xyKKf xh yhK 對對，即即四四個個點點，可可導(dǎo)導(dǎo)出出四四階階公公式式。常常用用的的四四 階階公公式式為為： 工程數(shù)學(xué)工程數(shù)學(xué)工程數(shù)學(xué)工程數(shù)學(xué)29 h=0.2,x=0 x=12(01);(0)1.xyyxyy 設(shè)設(shè)取取步步長長從從直直到到用用四四階階龍龍格格庫庫塔塔方方法法求求解解初初值值問問題題例例：工程數(shù)學(xué)工程數(shù)學(xué)工程數(shù)學(xué)工程數(shù)學(xué)30112341211322433(22);62;2;222;222().nnnnnnnnnnnnnnhyy

17、KKKKxKyyxhhKyKhyKxhhKyKhyKxhKyh Kyh K 由由 經(jīng)經(jīng) 典典 的的 四四 階階 龍龍 格格 庫庫 塔塔解解公公 式式 得得：工程數(shù)學(xué)工程數(shù)學(xué)工程數(shù)學(xué)工程數(shù)學(xué)312RK RKRK RK）方方法法的的導(dǎo)導(dǎo)出出基基于于T Ta ay yl lo or r展展開開，故故要要求求所所求求問問題題的的解解具具有有較較高高的的光光滑滑度度。當當解解充充分分光光滑滑時時，四四階階方方法法確確實實優(yōu)優(yōu)于于改改進進E Eu ul le er r法法。對對一一般般實實際際問問題題，四四階階方方法法一一般般可可達達到到精精度度要要求求。如如果果解解的的光光滑滑性性差差，則則用用四四階階

18、方方法法解解的的效效果果不不如如改改進進E Eu ul le er r法法。兩點說明兩點說明:1RKRK46RK5）當當p=1,2,3,4p=1,2,3,4時時，公公式式的的最最高高階階數(shù)數(shù)恰恰好好是是p,p,當當p4p4時時，公公式式的的最最高高階階數(shù)數(shù)不不是是p p，如如p=5p=5時時仍仍為為 ，p=p= 時時公公式式的的最最高高階階數(shù)數(shù)為為 。工程數(shù)學(xué)工程數(shù)學(xué)工程數(shù)學(xué)工程數(shù)學(xué)32五、變步長的龍格五、變步長的龍格庫塔方法庫塔方法()1()51115(2)15(2)11(2)11()11,(),2,2()2,2()1.()16hnhnnnnhnhnnhnnhnnhyy xychhxxhyc

19、hy xycy xyy xy n n以以經(jīng)經(jīng)典典四四階階龍龍格格庫庫塔塔公公式式為為例例。從從節(jié)節(jié)點點x x 出出發(fā)發(fā)，以以 為為步步長長求求一一近近似似值值將將步步長長折折半半，即即取取為為步步長長從從跨跨兩兩步步到到，求求一一近近似似值值每每跨跨一一步步的的截截斷斷誤誤差差是是因因此此有有由由上上兩兩式式 (2)(2)()11111().15hhhnnnny xyyy 工程數(shù)學(xué)工程數(shù)學(xué)工程數(shù)學(xué)工程數(shù)學(xué)33R-K方法的絕對穩(wěn)定區(qū)域/2121233223443 11 ,22111 )22411 )24 ( , ) () K()()()()()nnnnnnnRKhhhhf x yyhhyyyyh

20、KKyyhhKKhyyhhhKK 代入公式：代入公式：將將234112341 2)6111 12624(2()()()nnnhyyKKKKyhhh 工程數(shù)學(xué)工程數(shù)學(xué)工程數(shù)學(xué)工程數(shù)學(xué)34234111112624()()()nnhhhh 則則234111 112624()()()hhhh 絕絕對對穩(wěn)穩(wěn)定定區(qū)區(qū)域域： 2 1 -3 -2 -1 0 -1 -2 工程數(shù)學(xué)工程數(shù)學(xué)工程數(shù)學(xué)工程數(shù)學(xué)3504年：設(shè)年：設(shè)f(x)四階連續(xù)可導(dǎo)，四階連續(xù)可導(dǎo)，試建立如下數(shù)值微分公式試建立如下數(shù)值微分公式 并推導(dǎo)該公式的截斷誤差。并推導(dǎo)該公式的截斷誤差。0,0,1,2.ixxih i01212( ) 2 ( )( )( )f xf xf xf xh234(4)011111234(4)21111222(4)012128.(10) ()( )( )( )( )( )23!4!()( )( )( )( )()23!4!()2 ( )()( )( )12hhhf xf xhfxfxfxfhhhf xf xhfxfxfxfhf xf xf xhfxfh分兩式相加除 ，得解：解：工程數(shù)學(xué)工程數(shù)學(xué)工程數(shù)學(xué)工程數(shù)學(xué)3605年：設(shè)年：設(shè)f(x)四階連