計(jì)算機(jī)數(shù)值方法教學(xué)課件第二章 常微分方程數(shù)值解法

《計(jì)算機(jī)數(shù)值方法教學(xué)課件第二章 常微分方程數(shù)值解法》由會(huì)員分享，可在線閱讀，更多相關(guān)《計(jì)算機(jī)數(shù)值方法教學(xué)課件第二章 常微分方程數(shù)值解法（81頁(yè)珍藏版）》請(qǐng)?jiān)谘b配圖網(wǎng)上搜索。

1、第二章第二章 常微分方程數(shù)值解法常微分方程數(shù)值解法Chapter 2 Numerical Solution of Ordinary Differential Equation(s)22.1 引引 言言 Fm假設(shè)牽引力F為恒定值dvFdtm為了確定待定常數(shù)，可以給定初始條件：或者給定邊界條件：00tvvendtendvv,dvmF t vdt假設(shè)牽引力不恒定呢？求速度)(tv3雖然求解微分方程有許多解析方法，但解析方法只能夠求解一些特殊類型的方程。還有一類近似方法稱為數(shù)值方法，它可以給出解在一些離散點(diǎn)上的近似值。利用計(jì)算機(jī)解微分方程主要使用數(shù)值方法。解析方法與數(shù)值方法解析方法與數(shù)值方法4主要研究

2、對(duì)象：初值問(wèn)題 00)()(),(yxybxayxfdxdy00()(,()xxy xyf x y xdx求數(shù)值解：求數(shù)值解：求求y(x)在離散數(shù)據(jù)點(diǎn)在離散數(shù)據(jù)點(diǎn)xk處的近似值處的近似值yk。y=y(x)xyx0=ax1x2x3xk-1xn-1xn=b xk5則稱 f(x,y)對(duì)y 滿足李普希茲條件，此時(shí)初值問(wèn)題在a,b上存在唯一的連續(xù)可微的解唯一的連續(xù)可微的解。定理1：設(shè) f(x,y)是定義在區(qū)域 G=(x,y)|axb,yR上的連續(xù)函數(shù)，若存在正的常數(shù) L 使：1212|(,)(,)|f x yf x yL yy12,xa byy使得對(duì)任意的及都成立，(Lipschitz)條件條件有解條件

3、：有解條件：6則Lipschitz條件成立：在 f(x,y)對(duì)y可微的情況下，若偏導(dǎo)數(shù)有界：(,),(,)f x yLx yGy有解條件的判斷：有解條件的判斷：*121212(,)|(,)(,)|()|f x yf x yf x yyyL yyy7定理2：如果f(x,y)在G=(x,y)|axb,yR上滿足Lipschitz條件，則初值問(wèn)題是適定的。適定性：指初值問(wèn)題中，初始值y0及微分方程的右端函數(shù)f(x,y)有微小變化時(shí)，只能引起解的微小變化。適定性條件：適定性條件：8在解的存在區(qū)間 a,b上取n+1個(gè)節(jié)點(diǎn) bxxxxan210這里把 1,.iiihxx i=0,1,n稱為由xi到xi+1

4、的步長(zhǎng) 一般取成等間距的：nabh求解方法：求解方法：步進(jìn)法步進(jìn)法（分為（分為單步法單步法和和多步法多步法）數(shù)值方法的基本思想數(shù)值方法的基本思想9本章規(guī)定：本章規(guī)定：在 處初值問(wèn)題的理論解用 表示，數(shù)值解法的近似解用 表示。nx()ny xny記 ，它和 是不同的，后者等于 。(,)nnnff xy(,()nnf xy x()ny x102.2 幾種簡(jiǎn)單的數(shù)值方法幾種簡(jiǎn)單的數(shù)值方法（一）（一）歐拉（歐拉（Euler）法）法 00)()(),(yxybxayxfdxdy001()(,)0,1,2,.nnnnyy xyyhf xyn 11、泰勒公式解釋、泰勒公式解釋、求導(dǎo)的兩點(diǎn)公式解釋、求導(dǎo)的兩點(diǎn)

5、公式解釋、積分公式解釋、積分公式解釋歐拉公式的的分析解釋歐拉公式的的分析解釋001()(,)nnnny xyyyhf xy 00(,)()yf x yy xy 121)1()(21)!1()(!)(!2)()()()(pppnpnnnnhpyhPxyhxyhxyxyxy2()(,()()nnny xhf xy xO h泰勒公式解釋泰勒公式解釋其中：1nnhxx可以得到：001()(,)nnnny xyyyhf xy 13求導(dǎo)的兩點(diǎn)公式解釋求導(dǎo)的兩點(diǎn)公式解釋1()()()()(,)nnnnnny xy xyxhyxf xy可以得到：可以得到：001()(,)nnnny xyyyhf xy 140

6、0(,)(1.1)()(1.2)yfx yaxby xy 對(duì)微分方程對(duì)微分方程(1.1)(1.1)兩端從兩端從1nnxx到進(jìn)行積分進(jìn)行積分11(,)nnnnxxxxy dxfx y dx11()()(,)nnxnnxy xy xf x y dx積分公式解釋積分公式解釋1511)(,()(,()nnnnnnnnxxf xy xyyf xy xh右端積分用左矩形數(shù)值求積公式：右端積分用左矩形數(shù)值求積公式：11()()(,)nnxnnxy xy xf x y dx001()(,)nnnny xyyyhf xy 即：即：(,)f x yfxnx1nx16歐拉公式的的幾何描述yxx0 x1 x2 x3

7、x4y=y(x),.2,1,0),()(100nyxhfyyyxynnnn00(,)()yf x yy xy 00,xy 11,xy17例題例題1：（取步長(zhǎng)：（取步長(zhǎng)h=0.1）用）用Euler方法求滿足條方法求滿足條件件 的的y(t)數(shù)值解。數(shù)值解。22,12(1.0)0.0tdyyt etdtty )2(21ntnnnnnetythyy 解：解：2718281828.0)2(1.00200001 tetytyy122111120.1()0.684755578tyyyt et 18ntnyny(tn)y(tn)-yn01234101.01.11.21.31.43.00.00.271830.6

8、84761.276983.0935515.398240.00.345920.866641.607223.6203618.6830.00.074090.181880.330240.526813.28486數(shù)值解列表為數(shù)值解列表為)(2eetyt 22,12(1.0)0.0tdyyt etdtty )2(21ntnnnnnetythyy 19歐拉方法的誤差估計(jì)歐拉方法的誤差估計(jì)()-nnney xy 通過(guò)數(shù)值方法進(jìn)行計(jì)算時(shí)，考慮每一步產(chǎn)生的誤差，通過(guò)數(shù)值方法進(jìn)行計(jì)算時(shí)，考慮每一步產(chǎn)生的誤差，從從x0開(kāi)始一步步累積到開(kāi)始一步步累積到xn，稱，稱 為該數(shù)值方法在為該數(shù)值方法在xn點(diǎn)處的點(diǎn)處的整體截?cái)嗾`

9、差整體截?cái)嗾`差，該誤差，該誤差與與xn 及之前的各步計(jì)算誤差都有關(guān)系。及之前的各步計(jì)算誤差都有關(guān)系。20歐拉方法的誤差估計(jì)歐拉方法的誤差估計(jì)*1*1111()(,(),)()-()-()(,(),)nnnnnnnnnnnyy xhxy xhRy xyy xy xhxy xh 為了簡(jiǎn)化分析，著重分析為了簡(jiǎn)化分析，著重分析xn點(diǎn)單步計(jì)算產(chǎn)生的誤差，點(diǎn)單步計(jì)算產(chǎn)生的誤差，即把即把xn點(diǎn)之前的計(jì)算當(dāng)作無(wú)誤差：點(diǎn)之前的計(jì)算當(dāng)作無(wú)誤差：稱該誤差為數(shù)值方法在稱該誤差為數(shù)值方法在xn+1點(diǎn)處的點(diǎn)處的局部截?cái)嗾`差局部截?cái)嗾`差。12111(,()()ppnnnRH xy xhO h局部截?cái)嗾`差的第一個(gè)非零項(xiàng)為局部

10、截?cái)嗾`差的第一個(gè)非零項(xiàng)為局部截?cái)嗾`差主項(xiàng)局部截?cái)嗾`差主項(xiàng)。21歐拉方法的誤差估計(jì)歐拉方法的誤差估計(jì) 如果求解公式的局部截?cái)嗾`差為如果求解公式的局部截?cái)嗾`差為R(h)=O(hp+1)，則則稱該求解公式稱該求解公式具有具有p階精度階精度，稱該方法為稱該方法為p階方法階方法。定義：定義：歐拉方法：歐拉方法：212()()()()2!()(,()()nnnnnnny xy xy xy x hhy xhf xy xO h*1()(,()nnnnyy xhf xy x*2111()-()nnnRy xyO h具有具有1階精度。階精度。22（二）向后歐拉法（二）向后歐拉法 hxyxyxynnn)()()(1

11、1 11100(,),0,1,2,(),nnnnyyhf xyny xy （1）方法）方法 其公式為：其公式為：()(,)y xf x yxfnx1nx23（2）局部截?cái)嗾`差）局部截?cái)嗾`差*232111()()()()2nnnnyxRy xyhO hO h 231()()()()()2nnnny xy xy xy x hhO h*11()()nnnyy xy xh11100(,),0,1,2,(),nnnnyyhf xyny xy 21()()()()nnny xy xy x hO h24例題例題2：用向后：用向后Euler法解初值問(wèn)題法解初值問(wèn)題 1)0()10(,2yxyxyy11112(

12、)nnnnnxyyh yy 向后向后Euler法的公式為法的公式為 解：解：x0=0,y0=1,取h=0.125方法比較及推廣方法比較及推廣：nEuler方法方法 顯式公式顯式公式n向后向后Euler方法方法 隱式公式隱式公式 n解一個(gè)非線性方程解一個(gè)非線性方程 難求解難求解 n顯式和隱式相結(jié)合顯式和隱式相結(jié)合 隱式的顯化隱式的顯化 01110()(,),nnnnxyyyhfy xy 100(),)nnnnxyyyhfyy x 26計(jì)算公式為：計(jì)算公式為：001111)(,.2,1,0),(),(yxynyxhfyyyxhfyynnnnnnnn1 ny 由顯式得到，稱為由顯式得到，稱為預(yù)估值；

13、預(yù)估值；yn+1由隱式得到，稱為由隱式得到，稱為校正值。校正值。這種求解方法統(tǒng)稱為這種求解方法統(tǒng)稱為預(yù)估校正方法預(yù)估校正方法。其求解過(guò)。其求解過(guò)程為：程為：nnyyyyyyy2211027例3 用預(yù)估校正方法求解微分方程（取h=0.1）：1)0()10(,2yxyxyy解：1000.1)2(00001 yxyhyy0918.1)1000.11.021000.1(1.01)2(11101 yxyhyy1827.1)0918.11.020918.1(1.00918.1)2(11112 yxyhyy1763.1)1827.12.021827.1(1.00918.1)2(22212 yxyhyy232

14、222()1.2599xyyh yy 332332()1.2547xyyh yy 28（三）梯形公式（三）梯形公式 00)()(),(yxybxayxfdxdy 11,nnxnnxy xy xfx y xdx 111,2nnnnnnhyyfxyfxy()(,)y xf x yfxnx1nx29梯形公式局部截?cái)嗾`差梯形公式局部截?cái)嗾`差 *1111(),()()()22nnnnnnnnnhhyy xfxyfxyy xy xy x231()()()()()2!3!nnnnny xyxy xy xy x hhh21()()()()2!nnnnyxy xy xy x hh*111()nnnRy xy 3

15、3112nyxhh 30 111,2nnnnnnhyyfxyfxy預(yù)估-校正方法：稱為改進(jìn)的Euler求解公式或改進(jìn)Euler法。11,2nnnnnnnnhyyfxyfxyhfxy111100(,)(,)(,),0,1,2,2()nnnnnnnnnnyyhf x yhyyf x yf xyny xy 31為了表示方便，可以改寫為：112121h()2(,)(,)nnnnnnyyKKKf xyKf xh yhK 11,2nnnnnnnnhyyfxyfxyhfxy32（四）歐拉方法的收斂性分析（四）歐拉方法的收斂性分析由初值問(wèn)題的單步法產(chǎn)生的近似解 ，如果對(duì)于任一固定的 均有 ，則稱該方法是收斂的

16、。0nxxnh0lim()nnhnyy x ny定義：局部截?cái)嗾`差：局部截?cái)嗾`差：若初值問(wèn)題的一個(gè)單步法的局部截?cái)嗾`差為：定理：11(),1,pnRO hp 且增量函數(shù)(,(),)nnxy xh 關(guān)于y滿足Lipschitz條件，則整體截?cái)嗾`差：111()-()pnnney xyO h整體截?cái)嗾`差比局整體截?cái)嗾`差比局部截?cái)嗾`差低部截?cái)嗾`差低1 1階階*1111()-()-()(,(),)nnnnnnnRy xyy xy xhxy xh 33證明：存在常數(shù)c，使得*111()Pnny xych*1111()()PnnnRy xyO h*11()(,(),)(,)(1)()nnnnnnnnnnyy

17、y xyhxy xhxyhhL y xy*1()(,(),)nnnnyy xhxy xh 1(,)nnnnyyhxyh *111111111111211210111()()(1)(1)(1)1(1)(1)1(1)(1).(1)(1)(1)1(1)(1)1nnnnnnnPnPPnPnPnnnPney xyy xyyychhL echhL chhL echhLhLechhLhLhLhLehLchhLehL 03411110(1)1(1)(1)1nPnnhLechhLehL 000(),0y xye21011().2hLhLhLhLe1(1)0(1)nnhLhLe(1)11PnhLnceheL 當(dāng)

18、固定時(shí)，1nxx 10(1)nnhxxba 所以()111Pb a Lpncehec hL 35（五）歐拉方法的穩(wěn)定性分析（五）歐拉方法的穩(wěn)定性分析問(wèn)題：常微分方程初值問(wèn)題數(shù)值解的每步計(jì)算都是在前一步計(jì)算的結(jié)果上進(jìn)行的，所以必須考慮前面的誤差對(duì)以后計(jì)算結(jié)果的影響，誤差的積累會(huì)不會(huì)蓋過(guò)真解呢？選用代表性試驗(yàn)方程：y=y (Re()0，在計(jì)算yn時(shí)引入了誤差n。若這個(gè)誤差在計(jì)算后面的yn+k(k=1,2,.)中所引的誤差n+k按絕對(duì)值均不增加，就說(shuō)這個(gè)數(shù)值方法對(duì)于這個(gè)步長(zhǎng)h和復(fù)數(shù) 是絕對(duì)穩(wěn)定的。若在區(qū)域R內(nèi)數(shù)值方法是絕對(duì)穩(wěn)定的，則R為該數(shù)值方法的絕對(duì)穩(wěn)定區(qū)域(區(qū)間)。36把歐拉方法用于試驗(yàn)方程：y

19、=y1nnnyyhy 1nnnh 誤差方程：11nnh 11h 要求誤差不增加：O-2-1Re(h)Im(h)37把向后歐拉方法用于試驗(yàn)方程：y=y111nnh O21Re(h)Im(h)11h 要求誤差不增加：可見(jiàn)隱式的向后Euler方法比顯式的Euler方法的絕對(duì)穩(wěn)定域要大得多。同階精度的數(shù)值方法，往往隱式方法比顯式方法的絕對(duì)穩(wěn)定域大。38例4：以 y=y為例判斷梯形公式的穩(wěn)定性：11()2nnnnhyyyy解出：11/21/2nnhyyh 1/211/2hh 0 這種穩(wěn)定性稱為無(wú)條件穩(wěn)定！39例5：用Euler法、向后法、向后Euler法、法、改進(jìn)的歐拉法（梯形公式）解初值問(wèn)題 83,(

20、12)(1)3yyxy 取步長(zhǎng)h=0.2，小數(shù)點(diǎn)后至少保留4位。解：Euler法 1831.60.4nnnnyyhyy 向后Euler法 1180.625113nnnyyhyh 1183nnnyyhy40改進(jìn)的Euler法13161371 nnyy111(,)(,)2nnnnnnhyyf xyf xy 383822.011 nnnnyyyy梯形公式法 11,2nnnnnnnnhyyfxyfxyhfxy 11624690.581.122nnnnnhyyhyhyy 410013xyh=0.23.00003.30083.46593.55653.60623.63351.01.21.41.61.83.0

21、真解真解改進(jìn)的歐拉法梯形公式梯形公式向后向后EulerEulerxi 3.0000 3.0000 3.0000 3.0000 3.4000 3.2500 3.3077 3.2800 3.5600 3.4063 3.4734 3.4424 3.6240 3.5039 3.5626 3.5366 3.6496 3.5649 3.6106 3.5912 3.6598 3.6031 3.6365 3.622913161371 nnyy11.60.4nnyy 10.6251nnyy 10.581.12nnyy 如何得到高精度的求解公式？422.3 Runge-Kutta法法 如何構(gòu)造更高精度的求解公式？

22、如何構(gòu)造更高精度的求解公式？一種思路是：對(duì)微分方程右端積分采用高次多項(xiàng)式近似，如用二次多項(xiàng)式近似可得到 )(,()(,(4),(62)(,()()(112222nnnnnnxxnnxyxfxyxfyxfhdxxyxfxyxynn 43 ),(),(4),(311222nnnnnnnnyxfyxfyxfhyy 該公式稱為Simpson公式。該公式的局部截?cái)嗾`差為：4(4)4()()90khRyO h 44 rnrnnnnyrhyhyhyy!221 Taylor展開(kāi)：(,)nnyf xyf (,)(,)xnnynnnxyyfxyfxyyfff 222xxxyxyyyyfffffffyff 另一種得

23、到高階方法的想法是直接利用泰勒級(jí)數(shù)展開(kāi)。如果能計(jì)算得到y(tǒng)的高階微商,則可寫出r階的計(jì)算公式45例題：取h=0.1，求解初值問(wèn)題 2,(0)1;00.5yyyx 解：2324(4)35,22,66,2424yyyy yyyyyyyyyy2231(2)2nnnnhyyhyy 23423451(2)(6)(24)2!3!4!nnnnnnhhhyyhyyyy 一階公式：二階公式：四階公式：21(1)nnnnnyyhyyhy 461()1y xx 真解：p0.10.20.30.40.511.100001.221001.370081.557791.8004621.110001.246891.421741.

24、662621.9208741.111101.249661.428481.666451.99942y(xn)1.111111.250021.428571.666672.00000 xn問(wèn)題：1、求微商麻煩；2、計(jì)算量大。47R-K方法是通過(guò)對(duì)不同點(diǎn)上的函數(shù)值做線性組合，構(gòu)造近似公式，把近似公式和泰勒展開(kāi)相比較，使前面的若干項(xiàng)相吻合，從而使近似公式達(dá)到一定的階數(shù)。問(wèn)題：如何組合函數(shù)值？NiiinnKcyy11),(1nnyxhfK -11(,)iininijjjKhf xa h yb K其中，48選擇參數(shù)ci,ai,bij 的原則是，要求的Runge-Kutta法的右端項(xiàng)在(xn,yn)處泰勒展開(kāi)

25、后按h的冪次重新整理得到的結(jié)果231123nnyyd hd hd h 與微分方程的解y(xn+1)在xn處的Taylor展開(kāi)式 2311.2!3!nnnnny xhy xhyxh yxh yx 有盡可能多的項(xiàng)重合。49以計(jì)算兩個(gè)函數(shù)值為例說(shuō)明 11122122211,nnnnnnyyc hKc hKKfxyKfxa h yb hK 選擇c1、c2、a2、b21*111R()nnny xy階最高 1,nnnKfxy xyx 2221122211,nnnnxyKfxa h y xb hKfxy xa hfb hK fO h*1ny 為了計(jì)算50*1122321222()nnnxyyy xcchyx

26、bc a hff fO ha于是 )(!21321hOxyhxyhxyxynnnn 若要求局部截?cái)嗾`差達(dá)到)(3hO,則要求有 12222121,1/2,/1ccc aba 122211,12ccab 選取51122211,12ccab 1121211122,nnnnnnyyhKhKKfxyKfxh yhK 二階R-K方法的公式 它的局部階段誤差為O(h3)。這是計(jì)算兩次函數(shù)值的情況下所能達(dá)到的最高階。11122122211,nnnnnnyyc hKc hKKfxyKfxa h yb hK 521222101,1/2ccab 取取，12121,11,22nnnnnnyyhKKfxyKfxh y

27、hK 中間點(diǎn)法：12221132,443ccab 取取，1121211344,22,33nnnnnnyyhKhKKfxyKfxh yhK 二階休恩(Heun)法：53經(jīng)典的R-K方法是一個(gè)四階的方法,公式為 112341213243226,/2,/2/2,/2,nnnnnnnnnnhyyKKKKKfxyKfxhyhKKfxhyhKKfxh yhK 1、一步法，可以自開(kāi)始；特點(diǎn)：2、精度較高；3、便于改變步長(zhǎng)；4、計(jì)算量較小。54例2：用二階R-K方法和四階R-K方法(h=0.1)求解解：,101fx yxy000,0,0.1,0.5axyhb 1121212101101nnnnnnhyyKKK

28、xyKxhyhK 二階R-K方法 0)0()10()1(10yxyxy 11234121324322/610110/21/210/21/2101nnnnnnnnnnyyh KKKKKxyKxhyhKKxhyhKKxhyhK 四階R-K方法55xn0.10.20.30.40.50.60.70.80.91.0二階二階R-K 0.0500 0.0012 0.0488 0.0000 0.0488 0.1830 0.0017 0.1813 0.0000 0.1813 0.3627 0.0004 0.3624 0.0000 0.3624 0.5475 -0.0031 0.5507 -0.0000 0.55

29、07 0.7059 -0.0076 0.7134 -0.0001 0.7135 0.8235 -0.0112 0.8346 -0.0001 0.8347 0.9012 -0.0125 0.9135 -0.0002 0.9137 0.9476 -0.0116 0.9590 -0.0002 0.9592 0.9733 -0.0093 0.9823 -0.0002 0.9826 0.9866 -0.0066 0.9931 -0.0002 0.9933誤差四階四階R-K誤差y=1-e-5x2說(shuō)明：1、高階方法比低階方法精度高，但需要y(x)光滑性好；2、步長(zhǎng)h并非越小越好；步長(zhǎng)選擇從大到小。56數(shù)學(xué)家

30、Butcher于1965年證明了運(yùn)算量與可以達(dá)到的最高精度階數(shù)關(guān)系如下：每步計(jì)算每步計(jì)算kR的個(gè)數(shù)的個(gè)數(shù)2345678910可以達(dá)到最可以達(dá)到最高精度階數(shù)高精度階數(shù)23445666R-2可見(jiàn)，通常只選用可見(jiàn)，通常只選用R 4的公式是合適的。的公式是合適的。57自動(dòng)選步長(zhǎng)方法自動(dòng)選步長(zhǎng)方法 設(shè)在xn處之前誤差不計(jì)：)(nnxyy 誤差：(/2)(/2)()11111()15hhhnnnny xyyy ()511()hnnny xyC h(/2)511()2(/2)hnnny xyCh(/2)(/2)()11111()()15hhhnnnny xyyy 58隱式隱式R-KR-K方法方法類似于顯式類

31、似于顯式R-K公式，稍加改變，就得到公式，稍加改變，就得到隱式隱式R-K方法方法。11Lnniiiyyhk1,Lininii jjjkfxc h yc ha k顯式公式中對(duì)系數(shù)求和的上限是顯式公式中對(duì)系數(shù)求和的上限是i-1。而在隱式公式中。而在隱式公式中對(duì)系數(shù)求和的上限是對(duì)系數(shù)求和的上限是L，需要用迭代法求出，需要用迭代法求出Ki。推導(dǎo)。推導(dǎo)隱式公式的思路和方法與顯式隱式公式的思路和方法與顯式R-KR-K法類似。通常，同法類似。通常，同級(jí)的隱式公式可以獲得比顯式公式更高的階。級(jí)的隱式公式可以獲得比顯式公式更高的階。59通常，同級(jí)的隱式公式獲得比顯式公式更高的階。常通常，同級(jí)的隱式公式獲得比顯式

32、公式更高的階。常用的隱式用的隱式R-K法有法有:2,2111nnnnnyyhxhfyy1級(jí)級(jí)2階中點(diǎn)公式階中點(diǎn)公式:),(),(2111nnnnnnyxfyxfhyy2級(jí)級(jí)2階梯形公式：階梯形公式：n+1n121nn122nn12hy=y+(k+k)23+3h3+23k=fx+h,y+k+hk64123-33-23hk=fx+h,y+hk+k6124 2級(jí)級(jí)4階階R-K公式：公式：60隱式龍格隱式龍格-庫(kù)塔法庫(kù)塔法 ),.,1().,(.11111mjhKhKyhxfKKKhyymmjjijijmmii )2,2(1111KhyhxfKhKyyiiii其中其中2階方法階方法 的絕對(duì)穩(wěn)定區(qū)域?yàn)榈?/p>

33、絕對(duì)穩(wěn)定區(qū)域?yàn)?ReImg而而顯式顯式 1 4 階方法的絕對(duì)階方法的絕對(duì)穩(wěn)定區(qū)域?yàn)榉€(wěn)定區(qū)域?yàn)閗=1k=2k=3k=4-1-2-3-123ReImg無(wú)條件穩(wěn)定無(wú)條件穩(wěn)定龍格龍格-庫(kù)塔法穩(wěn)定區(qū)域庫(kù)塔法穩(wěn)定區(qū)域612.4 線性多步法線性多步法 單步法的優(yōu)缺點(diǎn)：簡(jiǎn)單，可以自開(kāi)始；提高精度時(shí)需要增加中間函數(shù)值；沒(méi)有充分利用前幾步得到的信息。多步法：多步法：yn-p,yn-p+1,yn-1 ,yn yn+162yn-p,yn-p+1,yn-1,yn yn+1的方法的方法考慮如下形式的求解公式11101kknin iin iiiyyhf此計(jì)算公式為線性的，所以稱為線性多步法。(,),(0,1,1)ijjif

34、f x yik(1,0,1)iik當(dāng) 時(shí)公式含有 ，這時(shí)公式是隱式的，而當(dāng) 時(shí)公式是顯式的011nf1063常微分初值問(wèn)題與積分公式等價(jià)11()()(,()nnxnnxy xy xf x y xdx 基本思想：基本思想：對(duì)f(x,y(x)進(jìn)行多項(xiàng)式插值，利用插值公式計(jì)算右邊積分，可以得到常微分初值問(wèn)題求解公式。64三次多項(xiàng)式插值求積分3121311231213112316121216()()()()()()()()()()()()()()()()()()nnnnnnnnnnnnnnnnF tLttxtxtxFxhtxtxtxFxhtxtxtxFxhtxtxtxFxh插值余項(xiàng)為插值余項(xiàng)為(4)3

35、11221()()()()()(),4!nnnnnnFRttxtxtxtxxx 令令 ，假設(shè)已知，假設(shè)已知 ，取插值節(jié)，取插值節(jié)點(diǎn)點(diǎn) ，對(duì)函數(shù)，對(duì)函數(shù) 作三次作三次Lagrange多項(xiàng)式插值多項(xiàng)式插值()(,()F xf x y x21,nnnyyy211,nnnnxxx x()F x651133()()()()nnxnnxy xy xLxRx dx 122110,nnxnjnnn in ixijn injj ixxyyfxydxxx 四階Adams內(nèi)插公式)5199(242111 nnnnnnffffhyy內(nèi)插公式局部截?cái)嗾`差：15(5)*32119(),720nnxnnxRx dxh yx

36、x 661133()()()()nnxnnxy xy xLxRx dx 令令 ，假設(shè)已知，假設(shè)已知 ，取插值節(jié)，取插值節(jié)點(diǎn)點(diǎn) ，對(duì)函數(shù)，對(duì)函數(shù) 作三次作三次Lagrange多項(xiàng)式插值多項(xiàng)式插值()(,()F xf x y x321,nnnnyyyy321,nnnnxxxx()F x1123(5559379)24nnnnnnhyyffff外插公式局部截?cái)嗾`差：15(5)331251(),720nnxnnxRx dxh yxx 133100,nnxnjnnn in ixijn injj ixxyyfxydxxx 67)5199(242111 nnnnnnffffhyy內(nèi)插公式1123(555937

37、9)24nnnnnnhyyffff外插公式通常計(jì)算公式為（預(yù)估校正系統(tǒng)）：)519),(9(24)9375955(242111113211nnnnnnnnnnnnnnfffyxfhyyffffhyy特點(diǎn)：精度高，但須與同精度的單步法配合使用。68例3：用Adams預(yù)測(cè)校正系統(tǒng)求初值問(wèn)題 3,00.0500dIEIIdtLLLtI 其中，E=200,L=3,=100,=50解：由于Adams預(yù)測(cè)校正系統(tǒng)為四階公式，選用四階R-K公式計(jì)算開(kāi)始值。a=t0=0;b=0.05I0=0;h=0.002f(t,I)=(200-50I3-100I)/3 121324311234,/2,/2/2,/2,22/

38、6nnnnnnnnnnKhf tIKhf thIKKhf thIKKhf th IKIIKKKK 691230.1289680.2493910.361332III再利用Adams預(yù)測(cè)校正系統(tǒng)1123111112(5559379)24(9(,)195)24nnnnnnnnnnnnnnhIIffffhIIfxIfff n45625t0.0080.0100.0120.050I0.4646910.5592060.644933.1.164861702.5 常微分方程組、高階微分方程及常微分方程組、高階微分方程及邊值問(wèn)題的數(shù)值解邊值問(wèn)題的數(shù)值解 考慮如下含二個(gè)未知函數(shù)的方程組：一、常微分方程組的數(shù)值解一、

39、常微分方程組的數(shù)值解 1112101022122020,(),()yfx yyyxyyfx yyyxy 00,yfx yy xy 則：可依照常微分方程數(shù)值方法來(lái)求解。令：1011120202212,yyfx yyyyfx yyyfx yy710000,yf x y zy xyzg x y zz xz 則其改進(jìn)的Euler格式具有預(yù)估形式11,nnnnnnnnnnyyhf xyzzzhg xyz校正公式為11111111,2,6nnnnnnnnnnnnnnnnhyyf xyzf xyzhzzg xyzg xyz如對(duì)于方程組如對(duì)于方程組72 1112111222121111222211221121

40、1321111222231122,/2,/2,/2/2,/2,/2/2,/2,/2/2,nnnnnnnnnnnnnnnnnnnKfxyyKfxyKfxyyKfxhyhKyhKKfxhyhKyhKKfxhyhKyhKKfxhyhK 24221111222211221,1133433311112,1122224322124/2,/2,226nnnnnnnnnnnyhKKfxh yhKyhKKfxh yhKyhKyyKKKKhyyKKKK 兩方程微分方程組的四階R-K公式73二、高階微分方程的數(shù)值解二、高階微分方程的數(shù)值解 對(duì)于高階微分方程，可以把它化成微分方程組后再進(jìn)行求解。如 0000,yfx

41、y yy xyyxy 令：112yyyy 則：1210010212200020 ,yyyxyyyfx yyyxyxy 使用四階R-K公式可得求解公式。74四階R-K公式 2112222111222221112222211111121122322334334222,/2/2,/2,/2/2/2,/2,/2,nnnnnnnnnnnnnnnnyKfxyyKyhKKfxhyhKyhKKyhKKfxhyhKyhKKyhKKfxh yhKyhKKy 31,1111112,1222222214314226nnnnyKKKKhyyKKKK 75三、邊值問(wèn)題的數(shù)值解三、邊值問(wèn)題的數(shù)值解 最簡(jiǎn)單的邊值問(wèn)題二階微分

42、方程及其定解條件:()(,(),(),(),()yxf x y xy xaxby ay b 它的解法可以用數(shù)值微分公式代替其導(dǎo)數(shù)值，將其變?yōu)榇鷶?shù)方程，然后求解，這種求解方法通常稱為差分法。222()()2()2()()kkkkkkky xhy xhyxO hhy xhy xy xhyxO hh 76 1111202(,)(1,2,1)2,kkkkkkknyyyyyf x yknhhy ayy by 其中，其中，h=(b-a)/n，xk=x0+kh，k=1,2,n-1，x0=a如果f是非線性函數(shù)，那么差分方程也是非線性的。差分方程為：上述差分方程為n-1個(gè)方程組成的n-1元代數(shù)方程組。解出y1,

43、y2,yn-1即為微分方程的計(jì)算解。77特別的，對(duì)二階線性常微分方程()()()yp x yq x yr x 其差分方程形式為 1111222kkkkkkkkkyyyyypq yrhh )1(.)1(.2112.211222122212121211212222222221212nhnnhnnnnhnnhhprhrhrhprhyyyyqhpphqhqhppqh整理得 78邊值問(wèn)題的定解條件分類：第一類邊界條件(),()y ay b第二類邊界條件(),()y ay b第三類邊界條件01010000()(),()()0,0,0y ay ay by b79例4 求解下面的邊值問(wèn)題，取h=0.2。213

44、630,100112xyxyyxxyyy 解：建立差分方程 211112110213()63212()kkkkkkkkyyyyyxxyxhhybyhay 將k=0代入(a)，聯(lián)立(b)得：102(2.520.5)2yyc 80將k=1，2，3，4代入(a)，聯(lián)立(c)得：012342.5220000.5212.07550.9623000.0679 012.03330.988900.0200120.91550.016900011.97671.7724yyyyy 解線性方程組得：012341.01321.01671.06931.21881.5130yyyyy 81小 結(jié)n求微分方程（組）數(shù)值解的單步法，多步法，但還有許多其它方法。n病態(tài)問(wèn)題（剛性問(wèn)題）求解比較困難。一般采用隱式方法進(jìn)行求解。n實(shí)際計(jì)算中，方法的選取非常重要。f(x,y)不太復(fù)雜時(shí)，采用單步法（四階R-K方法）；但當(dāng)y(x)的光滑性不好時(shí)，采用低階單步法（二階R-K方法）；f(x,y)比較復(fù)雜時(shí)，采用多步法。