數(shù)值計(jì)算方法-第9章-常微分方程數(shù)值解

1、二二0一一0年九月一日年九月一日復(fù)旦大學(xué)力學(xué)與工程科學(xué)系復(fù)旦大學(xué)力學(xué)與工程科學(xué)系飛行器設(shè)計(jì)與工程專業(yè)飛行器設(shè)計(jì)與工程專業(yè)主講人：楊永明主講人：楊永明計(jì)算方法計(jì)算方法講義講義1第9章 常微分方程 數(shù)值解法第第9章常微分方程數(shù)值解法章常微分方程數(shù)值解法復(fù)旦大學(xué)力學(xué)與工程科學(xué)系復(fù)旦大學(xué)力學(xué)與工程科學(xué)系9-2很多實(shí)際問(wèn)題往往歸結(jié)為一組微分方程很多實(shí)際問(wèn)題往往歸結(jié)為一組微分方程, 但它們的解析解很難得但它們的解析解很難得到到, 往往采用數(shù)值方法得到近似解往往采用數(shù)值方法得到近似解. 本章主要討論常微分方程本章主要討論常微分方程(組組)初值問(wèn)題的數(shù)值解法初值問(wèn)題的數(shù)值解法.l本章的目的本章的目的概述概述u

2、 問(wèn)題的提法問(wèn)題的提法l常微分方程初值問(wèn)題的提法常微分方程初值問(wèn)題的提法假設(shè)假設(shè) 連續(xù)且關(guān)于連續(xù)且關(guān)于y滿足滿足李譜希茲李譜希茲(Lipschitz)條件條件:( , )f x y1212( ,)( ,)f x yf x yL yyL為給定的常數(shù)為給定的常數(shù).根據(jù)常微分方程理論根據(jù)常微分方程理論, 上面的方程一定存在唯一的連續(xù)可微解上面的方程一定存在唯一的連續(xù)可微解.求方程的解求方程的解 y=y(x).0( , )( )dyf x yaxbdxy ay (1)第第9章常微分方程數(shù)值解法章常微分方程數(shù)值解法復(fù)旦大學(xué)力學(xué)與工程科學(xué)系復(fù)旦大學(xué)力學(xué)與工程科學(xué)系9-3l邊值問(wèn)題的適定性邊值問(wèn)題的適定性在

3、許多實(shí)際問(wèn)題中在許多實(shí)際問(wèn)題中, 是由觀測(cè)得到的是由觀測(cè)得到的, 存在一定的誤差存在一定的誤差. 如果它們的誤差是微小的如果它們的誤差是微小的, 那么能否保證解的誤差也是微小的那么能否保證解的誤差也是微小的.0( , ),f x yy定義定義: 假設(shè)初值問(wèn)題假設(shè)初值問(wèn)題(1)有唯一解有唯一解 . 如果存在正常數(shù)如果存在正常數(shù) , 使得對(duì)任何使得對(duì)任何 , 當(dāng)當(dāng)( )yy x1,K100zy以及以及( ),xaxb攝動(dòng)的初值問(wèn)題攝動(dòng)的初值問(wèn)題:0( , )( ) (),( )zf x zxaxbz az有唯一解有唯一解z(x), 并且滿足并且滿足: 1( )( )z xy xK則稱初值問(wèn)題則稱初

4、值問(wèn)題(1)是是適定的適定的. 定理定理: 若若 連續(xù)連續(xù), 且關(guān)于且關(guān)于y滿足滿足Lipschitz條件條件, 則初值問(wèn)題則初值問(wèn)題(1)是適定的是適定的.( , )f x y第第9章常微分方程數(shù)值解法章常微分方程數(shù)值解法復(fù)旦大學(xué)力學(xué)與工程科學(xué)系復(fù)旦大學(xué)力學(xué)與工程科學(xué)系9-4u 離散變量法及離散誤差離散變量法及離散誤差l離散變量法離散變量法在若干個(gè)離散點(diǎn)在若干個(gè)離散點(diǎn): a=x0 x1xn=b 上上, 求出函數(shù)求出函數(shù)y(x)的近似值的近似值yi (i=1,2,n).稱稱yi為原問(wèn)題的為原問(wèn)題的數(shù)值解數(shù)值解(即用它來(lái)近似即用它來(lái)近似y(xi)的值的值). hi = xi+1-xi 為為xi

5、到到 xi+1的步長(zhǎng)的步長(zhǎng). 一般取步長(zhǎng)為常數(shù)一般取步長(zhǎng)為常數(shù)h.要求數(shù)值解要求數(shù)值解, 必須對(duì)微分方程進(jìn)行離散必須對(duì)微分方程進(jìn)行離散. 常用方法有以下幾種常用方法有以下幾種:10(,)( )nnnnyyhf xyyy a1()()()nnny xy xy xhl用差商近似導(dǎo)數(shù)用差商近似導(dǎo)數(shù)如用向前差商代替導(dǎo)數(shù)如用向前差商代替導(dǎo)數(shù)(, ()nnf xy x用用y(xn)的近似值的近似值yn代入上式代入上式, 可得可得y(xn+1)的近似值的近似值yn+1:1(,)nnnnyyhf xy原方程的近似解原方程的近似解:0,1,2n 也可用其它格式的差商代替導(dǎo)數(shù)也可用其它格式的差商代替導(dǎo)數(shù).第第9章

6、常微分方程數(shù)值解法章常微分方程數(shù)值解法復(fù)旦大學(xué)力學(xué)與工程科學(xué)系復(fù)旦大學(xué)力學(xué)與工程科學(xué)系9-5l用數(shù)值積分方法用數(shù)值積分方法用左端點(diǎn)的矩形公式計(jì)算積分用左端點(diǎn)的矩形公式計(jì)算積分11()()( , )nnxnnxy xy xf x y dx( , )dyf x ydx兩端兩端積分積分nx1nx( )f x1(,)nnnnyyhf xy與差商近似的公式相同與差商近似的公式相同得得y(xn+1)的近似值的近似值yn+1用梯形公式計(jì)算積分用梯形公式計(jì)算積分nx1nx( )f x111 (,)(,)2nnnnnnhyyf xyf xy111()() (, ()(, ()2nnnnnnhy xy xf xy

7、 xf xy x311(), ()12nnnnh yxx用用yn近似近似y(xn), yn+1近似近似y(xn+1), 得得梯形方法梯形方法第第9章常微分方程數(shù)值解法章常微分方程數(shù)值解法復(fù)旦大學(xué)力學(xué)與工程科學(xué)系復(fù)旦大學(xué)力學(xué)與工程科學(xué)系9-6l用用Talor 展開近似的方法展開近似的方法用用y(xn)的近似值的近似值yn代入上式代入上式, 可得可得y(xn+1)的近似值的近似值yn+1:2()11()( )( )( )( )2!ppy xhy xhy xh yxh yxp1(, )nnnnyyhxyh上面三種方法都能導(dǎo)出常微分方程數(shù)值解的公式上面三種方法都能導(dǎo)出常微分方程數(shù)值解的公式, 其中其中

8、Talor 展展開方法還可以用來(lái)估計(jì)截?cái)嗾`差開方法還可以用來(lái)估計(jì)截?cái)嗾`差, 所以推導(dǎo)時(shí)都用該方法所以推導(dǎo)時(shí)都用該方法.與差商近似的公式相同與差商近似的公式相同1(1)1( ),( ,)(1)!pphyx xhp記為記為( , , )hx y h1(1)11( , , )( , )( , )( , )2!ppx y hf x yhfx yhfx yp令令x=xn, 則則1(1)11()()(, (), )( )(1)!ppnnnny xy xhxy xhhyp若取若取p=1, 得得1(,)nnnnyyhf xy局部截?cái)嗾`差或局部截?cái)嗾`差或簡(jiǎn)稱截?cái)嗾`差簡(jiǎn)稱截?cái)嗾`差第第9章常微分方程數(shù)值解法章常微分

9、方程數(shù)值解法復(fù)旦大學(xué)力學(xué)與工程科學(xué)系復(fù)旦大學(xué)力學(xué)與工程科學(xué)系9-7一一. 歐拉歐拉 (Euler) 方法方法u Euler 方法的各種形式方法的各種形式lEuler 方法方法稱為稱為 Euler 方法方法10(,)0,1,2( )nnnnyyhf xynyy a近似解是通過(guò)近似解是通過(guò)(x0,y0)的一條折的一條折線線, 每個(gè)折線段的方向與左端點(diǎn)每個(gè)折線段的方向與左端點(diǎn)處處f(x)的切線方向一致的切線方向一致. 故故Euler 方法又稱為方法又稱為Euler折線法折線法.0( , )( )dyf x yaxbdxy ay數(shù)值解數(shù)值解lEuler 方法的幾何解釋方法的幾何解釋Euler 方 法

10、的 幾 何 解 釋P0P1P2P3P4Pn( )yf xOx0 x1 x2 x3 x4 xnxyEuler 方法是顯式的方法是顯式的, 可直接遞可直接遞推求解推求解.第第9章常微分方程數(shù)值解法章常微分方程數(shù)值解法復(fù)旦大學(xué)力學(xué)與工程科學(xué)系復(fù)旦大學(xué)力學(xué)與工程科學(xué)系9-8l向后向后Euler 方法方法向后向后Euler方法方法這是隱式公式這是隱式公式, 一般用迭代法求解一般用迭代法求解:若用向后差商代替導(dǎo)數(shù)若用向后差商代替導(dǎo)數(shù), 即即:211()()nnnyyhy xO h11()()()( )nnny xy xy xO hh用用y(xn)的近似值的近似值yn代入上式代入上式, 可得可得y(xn+1

11、)的近似值的近似值yn+1:111(,)nnnnyyhf xy(0)1(1)( )111(,)(,)(0,1,2,)nnnnkknnnnyyhf xyyyhf xykl中心中心Euler 方法方法若用中心差商代替導(dǎo)數(shù)若用中心差商代替導(dǎo)數(shù), 即即:211()()()()2nnny xy xy xO hh中點(diǎn)中點(diǎn)Euler方法方法3112()()nnnyyhy xO h112(,)nnnnyyhf xy因?yàn)橐玫角懊鎯刹降慕Y(jié)果因?yàn)橐玫角懊鎯刹降慕Y(jié)果yn-1,yn, 故又稱為故又稱為Euler兩步公式兩步公式.第第9章常微分方程數(shù)值解法章常微分方程數(shù)值解法復(fù)旦大學(xué)力學(xué)與工程科學(xué)系復(fù)旦大學(xué)力學(xué)與工程

12、科學(xué)系9-9lEuler 方法的算法描述方法的算法描述1. 輸入輸入a, b, h, f(x, y), 初值初值y0;2. n0, xa, yy0;3. 輸出輸出n, (x, y);4. nn+1, yy+hf(x, y);5. xx+h; 若若x 0, 使得使得整體截?cái)嗾`差整體截?cái)嗾`差:1111111()()nnnnnnney xyy xyyy1 ()(, ()(,)nnnnnnnRy xh f xy xyhf xy1()(, ()(,)nnnnnnnRy xyh f xy xf xy21()2nnnh MehL y xy21(1)2nh MhL e由閉區(qū)間上連續(xù)函數(shù)由閉區(qū)間上連續(xù)函數(shù)的性質(zhì)

13、的性質(zhì), Rn+1在在a,b區(qū)間有界區(qū)間有界.第第9章常微分方程數(shù)值解法章常微分方程數(shù)值解法復(fù)旦大學(xué)力學(xué)與工程科學(xué)系復(fù)旦大學(xué)力學(xué)與工程科學(xué)系9-14lEuler 方法的整體截?cái)嗾`差方法的整體截?cái)嗾`差 (續(xù)續(xù))反復(fù)遞推得反復(fù)遞推得,211(1)2nneh MhL e221111(1)(1)22nneh MhLh MhL e21001(1)(1)2nknkh MhLhLe021(1)12(1)1nh MhLhL1(1)12nhMhLL因?yàn)樵谝驗(yàn)樵赼,b中求解中求解, 故有故有:(1)()nhba(1)banh1(1)(1)b anhhLhL利用不等式利用不等式:1,(0)xxex()b ahLL

14、b ahee(取取x=hL)()112L b anhMeeLEuler方法的整體截?cái)嗾`差與方法的整體截?cái)嗾`差與h同階同階. 當(dāng)當(dāng)h0時(shí)時(shí), eN0第第9章常微分方程數(shù)值解法章常微分方程數(shù)值解法復(fù)旦大學(xué)力學(xué)與工程科學(xué)系復(fù)旦大學(xué)力學(xué)與工程科學(xué)系9-15二二. 改進(jìn)的歐拉改進(jìn)的歐拉 (Euler) 方法方法u 梯形公式梯形公式l梯形公式的導(dǎo)出梯形公式的導(dǎo)出梯形公式梯形公式l梯形公式的誤差梯形公式的誤差取取a=xn, b=xn+1, 得梯形方法的局部誤差得梯形方法的局部誤差:用梯形積分公式在用梯形積分公式在xn,xn+1上求積分上求積分:1111()()( , ) (, ()(, ()2nnxnnnn

15、nnxhy xy xf x y dxf xy xf xy x( , )yf x y111 (,)(,)2nnnnnnhyyf xyf xy用近似值用近似值yn代替代替y(xn), yn+1代替代替y(xn+1) , 可得可得:由定理由定理7.2梯形積分公式的誤差梯形積分公式的誤差:31()( )( ) ( )( )( ),212bababaRff x dxf af bf ( , )a b第第9章常微分方程數(shù)值解法章常微分方程數(shù)值解法復(fù)旦大學(xué)力學(xué)與工程科學(xué)系復(fù)旦大學(xué)力學(xué)與工程科學(xué)系9-16l梯形公式的誤差梯形公式的誤差 (續(xù)續(xù))梯形公式是隱式格式梯形公式是隱式格式, 通常采用迭代方法求解通常采用

16、迭代方法求解. 1()nnxx31111()() (, ()(, ()( )212nnnnnnnhhRy xy xf xy xf xy xf 3()O h梯形公式是二階方法梯形公式是二階方法l梯形公式的迭代格式梯形公式的迭代格式(0)1(1)( )111(,) (,)(,) (0,1,2,)2nnnnkknnnnnnyyh f xyhyyf xyf xykl梯形公式的收斂性梯形公式的收斂性由于由于f(x,y)關(guān)于關(guān)于y滿足滿足Lipschitz條件條件, (1)( )( )(1)111111(,)(,)2kkkknnnnnnhyyf xyf xy( )(1)112kknnhL yy(1)(0)

17、112knnhLyy若若 , 則迭代收斂則迭代收斂. 但這樣做計(jì)算量太大但這樣做計(jì)算量太大.21hLL是是Lipschitz常數(shù)常數(shù).第第9章常微分方程數(shù)值解法章常微分方程數(shù)值解法復(fù)旦大學(xué)力學(xué)與工程科學(xué)系復(fù)旦大學(xué)力學(xué)與工程科學(xué)系9-17u改進(jìn)改進(jìn)Euler 法法l改進(jìn)改進(jìn)Euler法的思想法的思想稱為稱為預(yù)測(cè)預(yù)測(cè)-校正系統(tǒng)校正系統(tǒng), 或稱或稱改進(jìn)改進(jìn)Euler法法. 它屬于顯式單步法它屬于顯式單步法.先用先用Euler法求法求yn+1的近似值的近似值 ;再用梯形公式迭代一次再用梯形公式迭代一次, 得到得到y(tǒng)n+1.1nyl編程采用的公式編程采用的公式1111(,) (,)(,)2nnnnnnn

18、nnnyyh f xyhyyf xyf xy校正校正預(yù)測(cè)預(yù)測(cè)11(,)(,)0.5()pnnnqnnpnpqyyh f xyyyh f xyyyy第第9章常微分方程數(shù)值解法章常微分方程數(shù)值解法復(fù)旦大學(xué)力學(xué)與工程科學(xué)系復(fù)旦大學(xué)力學(xué)與工程科學(xué)系9-18l改進(jìn)改進(jìn)Euler 法的算法描述法的算法描述1.輸入輸入a,b, f(x,y), 區(qū)間等分?jǐn)?shù)區(qū)間等分?jǐn)?shù)N, 初值初值y0;3.for i=1 to N step 1 計(jì)算各點(diǎn)的函數(shù)值計(jì)算各點(diǎn)的函數(shù)值4.停機(jī)停機(jī).2.計(jì)算步長(zhǎng)計(jì)算步長(zhǎng) h=(b-a)/N; 置置n0, xa, yy0; 輸出輸出(x,y);3.1 計(jì)算計(jì)算 yp=y+hf(x,y);

19、3.2置置xx+h; 計(jì)算計(jì)算yq=y+hf(x,yp);3.3置置y0.5(yp+yq);3.4輸出輸出(x,y);End for il改進(jìn)改進(jìn)Euler法的誤差法的誤差可以證明可以證明, 改進(jìn)改進(jìn)Euler 法的局部截?cái)嗾`差為法的局部截?cái)嗾`差為O(h3), 故屬于二階方故屬于二階方法法.第第9章常微分方程數(shù)值解法章常微分方程數(shù)值解法復(fù)旦大學(xué)力學(xué)與工程科學(xué)系復(fù)旦大學(xué)力學(xué)與工程科學(xué)系9-19l改進(jìn)改進(jìn)Euler 法的算例法的算例從結(jié)果可以看出從結(jié)果可以看出: 精確度比精確度比Euler法明顯提高法明顯提高.誤誤差差越越來(lái)來(lái)越越大大取步長(zhǎng)取步長(zhǎng)h=0.1; 結(jié)果如下結(jié)果如下:22, 12,(1)

20、0 xyy xx exy nxnyny(xn)y(xn)-yn01.000011.10.3423780.3459200.00354221.20.8583150.8666430.00832831.31.5927501.6072150.01446541.42.5982982.6203600.02206151.53.9364443.9676660.03122261.65.6789075.7209620.04205471.77.9092097.9638730.05466481.810.72446710.7936250.06915891.914.23744214.3230820.085640102.01

21、8.57888218.6830970.104215第第9章常微分方程數(shù)值解法章常微分方程數(shù)值解法復(fù)旦大學(xué)力學(xué)與工程科學(xué)系復(fù)旦大學(xué)力學(xué)與工程科學(xué)系9-20三三. 龍格龍格-庫(kù)塔庫(kù)塔 (Runge-Kutta) 法法u Runge-Kutta 法的基本思想法的基本思想l如何提高精度如何提高精度局部截?cái)嗾`差由局部截?cái)嗾`差由Talor展開式的余項(xiàng)決定展開式的余項(xiàng)決定, 一般地一般地,2()1()()()2!ppnnnnnhhyyhy xyxyxp具有具有p 階精度階精度, 局部截?cái)嗾`差為局部截?cái)嗾`差為:1()pO h其中其中, ( )( , )( , ),xyyxfx yfx y y( )( , ),

22、y xf x y( , )f x y特例特例. p=1時(shí)時(shí), 即為即為Euler方法方法.1(,)nnnnyyhf xy第第9章常微分方程數(shù)值解法章常微分方程數(shù)值解法復(fù)旦大學(xué)力學(xué)與工程科學(xué)系復(fù)旦大學(xué)力學(xué)與工程科學(xué)系9-21l障礙及對(duì)策障礙及對(duì)策如如, Euler 法可表示為法可表示為:改進(jìn)改進(jìn)Euler 法可表示為法可表示為:111(,)nnnnyyh KKf xy其其中中，112121(22),(,)(,)nnnnnnyyhKKKf xyKf xh yhK其其中中，障礙障礙: 必須求必須求f(x,y) 的各階偏導(dǎo)數(shù)的各階偏導(dǎo)數(shù), 計(jì)算復(fù)雜且工作量大計(jì)算復(fù)雜且工作量大.對(duì)策對(duì)策: 用用f(x,

23、y)在某些點(diǎn)上函數(shù)值的線性組合在某些點(diǎn)上函數(shù)值的線性組合,來(lái)計(jì)算近似值來(lái)計(jì)算近似值yn+1, 并使它的并使它的Talor展開式與展開式與y(xn+1)的展開式前面若干項(xiàng)完全相同的展開式前面若干項(xiàng)完全相同, 從而達(dá)到一定階數(shù)的精度從而達(dá)到一定階數(shù)的精度這就是這就是Runge-Kutta法法(簡(jiǎn)稱簡(jiǎn)稱RK方法方法) 的基本思想的基本思想.只計(jì)算一點(diǎn)的函數(shù)值只計(jì)算一點(diǎn)的函數(shù)值 , 得到的是一階方法得到的是一階方法.(,)nnf xy須計(jì)算兩點(diǎn)處的函數(shù)值須計(jì)算兩點(diǎn)處的函數(shù)值 , 其其Talor展開的前三項(xiàng)與展開的前三項(xiàng)與y(xn+1)的的展開式相同展開式相同, 得到的是二階方法得到的是二階方法.( ,

24、 )f x y第第9章常微分方程數(shù)值解法章常微分方程數(shù)值解法復(fù)旦大學(xué)力學(xué)與工程科學(xué)系復(fù)旦大學(xué)力學(xué)與工程科學(xué)系9-22uRKRK方法的構(gòu)造方法的構(gòu)造l構(gòu)造方法構(gòu)造方法1=111=1(,)(,)pnniiinniininijjjyyhc KKf xyKf xa h yhb K其其中中，一般地一般地, 設(shè)設(shè)RK近似公式為近似公式為:其中其中, 為待定參數(shù)為待定參數(shù).,iijia bcl二階二階RK方法方法(2,3,)ip11122122211()(,)(,)nnnnnnyyhc Kc KKf xyKf xa h yhb K取取p=2, 近似公式為近似公式為:在在(xn,yn)處處Talor展開展開:

25、222211(,)(,)(,)()nnxnnynnKf xyfxya hfxyhb KO h1122221(,) (,)(,)(,) (,)()nnnnnnxnnynnnnyyhc f xyhcf xyfxya hfxyf xyhbO h112232221()(,)(,)(,)(,)()nnnnxnnynnnnyyccf xyhc a fxybfxyf xyhO h整理得整理得:第第9章常微分方程數(shù)值解法章常微分方程數(shù)值解法復(fù)旦大學(xué)力學(xué)與工程科學(xué)系復(fù)旦大學(xué)力學(xué)與工程科學(xué)系9-23l二階二階RK方法（續(xù)）方法（續(xù)）理論上可以證明理論上可以證明, 無(wú)論怎樣選取參數(shù)無(wú)論怎樣選取參數(shù)c1, c2, a

26、2, b21, 上面的公式不可上面的公式不可能有更高的精度能有更高的精度. 即即, 每步計(jì)算兩個(gè)函數(shù)值每步計(jì)算兩個(gè)函數(shù)值, 只能得出二階公式只能得出二階公式.112232221()(,)(,)(,)(,)()nnnnxnnynnnnyyccf xyhc a fxybfxyf xyhO h1()ny x在在xn處處Talor展開展開23()()()()2nnnhy xhy xyxO h( )( , )y xf x y( )( , )( , )( )xyyxfx yfx y y x231()(,)(,)(,) (,)()2nnnnxnnynnnnhy xyhf xyfxyfxyf xyO hh的各

27、階的各階冪次相同冪次相同122222111 21 2ccc ac b共共4個(gè)未知量個(gè)未知量, 有無(wú)窮多個(gè)解有無(wú)窮多個(gè)解.代入代入RK近似公式近似公式, 其局部截?cái)嗾`差其局部截?cái)嗾`差都是都是O(h3). 統(tǒng)稱為統(tǒng)稱為二階方法二階方法.21222121 (2)1bacacc 第第9章常微分方程數(shù)值解法章常微分方程數(shù)值解法復(fù)旦大學(xué)力學(xué)與工程科學(xué)系復(fù)旦大學(xué)力學(xué)與工程科學(xué)系9-24l幾個(gè)常用的二階幾個(gè)常用的二階RK公式公式122222111 21 2ccc ac b改進(jìn)改進(jìn)Euler公式公式取取122211 2,1ccab11122122211()(,)(,)nnnnnnyyhc Kc KKf xyKf

28、 xa h yhb K中點(diǎn)公式中點(diǎn)公式取取122210,1,1 2ccabHeun二階公式二階公式取取122211 4,3 4,2 3ccab112121(3)4(,)22(,)33nnnnnnhyyKKKf xyKf xh yhKHeun 公式公式12121(,)(,)22nnnnnnyyhKKf xyhhKf xyK中點(diǎn)公式中點(diǎn)公式112121()2(,)(,)nnnnnnhyyKKKf xyKf xh yhK改進(jìn)改進(jìn)Euler公式公式第第9章常微分方程數(shù)值解法章常微分方程數(shù)值解法復(fù)旦大學(xué)力學(xué)與工程科學(xué)系復(fù)旦大學(xué)力學(xué)與工程科學(xué)系9-25l幾個(gè)常用的三階幾個(gè)常用的三階RK公式公式Kutta三

29、階公式三階公式RK一般公式中一般公式中, 取取p=3可得三階可得三階RK公式公式. 其中包含其中包含c1, c2, c3, a2, a3, b21, b31, b32 共共6個(gè)系數(shù)個(gè)系數(shù), 也有無(wú)窮多個(gè)解也有無(wú)窮多個(gè)解. 常用的三階公式有常用的三階公式有: Heun三階公式三階公式1123121312(4)6(,)(,)22(,2)nnnnnnnnhyyKKKKf xyhhKf xyKKf xh yhKhK11312132(3)4(,)(,)3322(,)33nnnnnnnnhyyKKKf xyhhKf xyKhhKf xyK第第9章常微分方程數(shù)值解法章常微分方程數(shù)值解法復(fù)旦大學(xué)力學(xué)與工程科學(xué)

30、系復(fù)旦大學(xué)力學(xué)與工程科學(xué)系9-26l幾個(gè)常用的四階幾個(gè)常用的四階RK公式公式經(jīng)典經(jīng)典Kutta四階公式四階公式RK一般公式中一般公式中, 取取p=4可得四階可得四階RK公式公式. 常用的四階公式有常用的四階公式有: 112341213243(22)6(,)(,)22(,)22(,)nnnnnnnnnnhyyKKKKKf xyhhKf xyKhhKf xyKKf xh yhK它是實(shí)際應(yīng)用中最常用的單步法它是實(shí)際應(yīng)用中最常用的單步法.第第9章常微分方程數(shù)值解法章常微分方程數(shù)值解法復(fù)旦大學(xué)力學(xué)與工程科學(xué)系復(fù)旦大學(xué)力學(xué)與工程科學(xué)系9-27l幾個(gè)常用的四階幾個(gè)常用的四階RK公式（續(xù)）公式（續(xù)）Gill公

31、式公式11234121312423(22)(22)6(,)(,)22212(,(1)22222(,(1)22nnnnnnnnnnhyyKKKKKf xyhhKf xyKhKf xyhKhKKf xh yhKhKlRK方法中方法中 p 值與階數(shù)的關(guān)系值與階數(shù)的關(guān)系p=1,2,3,4時(shí)時(shí), 得到得到RK公式的最高階數(shù)分別為一、二、三、四階公式的最高階數(shù)分別為一、二、三、四階.但最高階數(shù)并不一定等于但最高階數(shù)并不一定等于p, 當(dāng)當(dāng)p=5,6,7,8,9 時(shí)時(shí), RK公式的最高階公式的最高階數(shù)分別為數(shù)分別為: 4, 5, 6, 6, 7.第第9章常微分方程數(shù)值解法章常微分方程數(shù)值解法復(fù)旦大學(xué)力學(xué)與工程

32、科學(xué)系復(fù)旦大學(xué)力學(xué)與工程科學(xué)系9-28l經(jīng)典四階經(jīng)典四階Kutta 方法的算法描述方法的算法描述1.輸入輸入a,b, f(x,y), 區(qū)間等分?jǐn)?shù)區(qū)間等分?jǐn)?shù)N, 初值初值y0;3.for i=1 to N step 1 計(jì)算各點(diǎn)的函數(shù)值計(jì)算各點(diǎn)的函數(shù)值4.停機(jī)停機(jī).2.計(jì)算步長(zhǎng)計(jì)算步長(zhǎng) h=(b-a)/N; 置置xa, yy0; 輸出輸出(x,y);3.1 計(jì)算計(jì)算 K1= f(x,y);3.2置置xx+0.5*h; 計(jì)算計(jì)算:K2= f(x,y+0.5*h*K1);K3= f(x,y+0.5*h*K2); 3.3置置x a+i*h; 計(jì)算計(jì)算: K4= f(x,y+h*K3);3.4置置y y

33、 + h*K1+2(K2+K3)+K4 / 6;3.5輸出輸出(x,y);End for i3.for i=1 to N step 1 計(jì)算各點(diǎn)的函數(shù)值計(jì)算各點(diǎn)的函數(shù)值第第9章常微分方程數(shù)值解法章常微分方程數(shù)值解法復(fù)旦大學(xué)力學(xué)與工程科學(xué)系復(fù)旦大學(xué)力學(xué)與工程科學(xué)系9-29l經(jīng)典四階經(jīng)典四階RK RK 公式算例公式算例從結(jié)果可以看出從結(jié)果可以看出: 精確度比精確度比改進(jìn)改進(jìn)Euler法明顯提高法明顯提高.誤誤差差越越來(lái)來(lái)越越大大取步長(zhǎng)取步長(zhǎng)h=0.1; 結(jié)果如下結(jié)果如下:22, 12,(1)0 xyy xx exy nxnyny( xn)y( xn)- yn01.000011.10.3459100

34、.3459209.5892e-621.20.8666220.8666432.0843e-531.31.6071811.6072153.3731e-541.42.6203112.6203604.8245e-551.53.9676023.9676666.4396e-561.65.7208795.7209628.2201e-571.77.9637727.9638731.0169e-481.810.79350210.7936251.2288e-491.914.32293614.3230821.4581e-4102.018.68292718.6830971.7051e-4第第9章常微分方程數(shù)值解法章常微

35、分方程數(shù)值解法復(fù)旦大學(xué)力學(xué)與工程科學(xué)系復(fù)旦大學(xué)力學(xué)與工程科學(xué)系9-30u變步長(zhǎng)的變步長(zhǎng)的RKRK方法方法l為何要變步長(zhǎng)為何要變步長(zhǎng)由于由于y(x)是不均勻的是不均勻的, 采用等不長(zhǎng)計(jì)算時(shí)采用等不長(zhǎng)計(jì)算時(shí), 有些點(diǎn)上精度較高有些點(diǎn)上精度較高, 而而有些點(diǎn)則精度較低有些點(diǎn)則精度較低. 是否能根據(jù)精度要求是否能根據(jù)精度要求, 自動(dòng)調(diào)整當(dāng)前計(jì)算步的自動(dòng)調(diào)整當(dāng)前計(jì)算步的步長(zhǎng)呢步長(zhǎng)呢?局部截?cái)嗾`差局部截?cái)嗾`差O(hp+1), 即即( )111()hpnny xychnyp階方法階方法,步長(zhǎng)步長(zhǎng)h( )1hnyc與與 有關(guān)有關(guān), 假設(shè)假設(shè) 內(nèi)變內(nèi)變化不大化不大, 則可看作常數(shù)則可看作常數(shù).(1)1( ),p

36、nnyxx(1)1( ),pnnyxx在在. (1)nyp階方法階方法,步長(zhǎng)步長(zhǎng)h/2(2)1 2hnyp階方法階方法,步長(zhǎng)步長(zhǎng)h/2(2)1hny截?cái)嗾`差截?cái)嗾`差O(h/2)p+1, 即即1( /2)11()22phnnhy xyc. (2)l從誤差的分析說(shuō)起從誤差的分析說(shuō)起如果分兩步推進(jìn)如果分兩步推進(jìn),第第9章常微分方程數(shù)值解法章常微分方程數(shù)值解法復(fù)旦大學(xué)力學(xué)與工程科學(xué)系復(fù)旦大學(xué)力學(xué)與工程科學(xué)系9-31l變步長(zhǎng)控制方法變步長(zhǎng)控制方法-1-1具體實(shí)施步驟具體實(shí)施步驟( )( /2)( )1111(21) ()22phphphnnnny xyyy( )( /2)( )11112()21phhh

37、nnnnpy xyyy可根據(jù)可根據(jù)的大小來(lái)選擇合適的步長(zhǎng)的大小來(lái)選擇合適的步長(zhǎng).(2)2p- -(1)( /2)( )111(21) ()20pphhnnny xyy1.以步長(zhǎng)以步長(zhǎng)h, 從從xn點(diǎn)出發(fā)計(jì)算點(diǎn)出發(fā)計(jì)算y(xn+h) 的近似值的近似值 ;( )1hny. (3)2.以步長(zhǎng)以步長(zhǎng)h/2, 從從xn點(diǎn)出發(fā)點(diǎn)出發(fā), 用兩步計(jì)算用兩步計(jì)算y(xn+h) 的近似值的近似值 ;( /2)1hny3.如果如果(3)式的式的滿足要求的精度滿足要求的精度, 即即 , 且兩者相差不多且兩者相差不多, 則仍以則仍以h為步長(zhǎng)繼續(xù)計(jì)算為步長(zhǎng)繼續(xù)計(jì)算; 4.如果如果 , 說(shuō)明步長(zhǎng)太小說(shuō)明步長(zhǎng)太小. 反復(fù)以反

38、復(fù)以2h為步長(zhǎng)計(jì)算為步長(zhǎng)計(jì)算, 直到直到 為止為止. 這時(shí)這時(shí), 前一次所用的步長(zhǎng)前一次所用的步長(zhǎng)h即為合適的步長(zhǎng)即為合適的步長(zhǎng), 繼續(xù)推進(jìn)繼續(xù)推進(jìn); 5.如果如果 , 說(shuō)明步長(zhǎng)太大說(shuō)明步長(zhǎng)太大. 反復(fù)以反復(fù)以h/2為步長(zhǎng)計(jì)算為步長(zhǎng)計(jì)算, 直到直到 為止為止. 此時(shí)所用的步長(zhǎng)此時(shí)所用的步長(zhǎng)h即為合適的步長(zhǎng)即為合適的步長(zhǎng), 繼續(xù)推進(jìn)繼續(xù)推進(jìn); 第第9章常微分方程數(shù)值解法章常微分方程數(shù)值解法復(fù)旦大學(xué)力學(xué)與工程科學(xué)系復(fù)旦大學(xué)力學(xué)與工程科學(xué)系9-32l變步長(zhǎng)控制方法的改進(jìn)變步長(zhǎng)控制方法的改進(jìn)具體實(shí)施步驟的改進(jìn)具體實(shí)施步驟的改進(jìn)(1)與與(2)式相除式相除( )11( /2)11()2()hpnnhnn

39、y xyy xy即步長(zhǎng)減半后即步長(zhǎng)減半后, 誤差降為原來(lái)的約誤差降為原來(lái)的約1/2p倍倍. 由此可推得由此可推得, 若步長(zhǎng)為若步長(zhǎng)為 , 則誤差降為原來(lái)的約則誤差降為原來(lái)的約 倍倍.2 ,(1,2,)khk 1 2kp13步與前相同步與前相同;4.如果如果 , 說(shuō)明步長(zhǎng)太小說(shuō)明步長(zhǎng)太小. 將步長(zhǎng)放大至將步長(zhǎng)放大至2kh, 則誤差放大至則誤差放大至 . 令令2kp取滿足條件的最大整數(shù)取滿足條件的最大整數(shù)k, 下一步的步長(zhǎng)取為下一步的步長(zhǎng)取為2kh, 繼續(xù)推進(jìn)繼續(xù)推進(jìn);22kpkp ln()ln 2kp5.如果如果 , 說(shuō)明步長(zhǎng)太大說(shuō)明步長(zhǎng)太大. 將步長(zhǎng)縮小至將步長(zhǎng)縮小至h/2k, 則誤差縮小至則

40、誤差縮小至 . 令令 2kp22kpkpln()ln 2kp取滿足條件的最小整數(shù)取滿足條件的最小整數(shù)k, 下一步的步長(zhǎng)取為下一步的步長(zhǎng)取為h/2k, 繼續(xù)推進(jìn)繼續(xù)推進(jìn);第第9章常微分方程數(shù)值解法章常微分方程數(shù)值解法復(fù)旦大學(xué)力學(xué)與工程科學(xué)系復(fù)旦大學(xué)力學(xué)與工程科學(xué)系9-331=0= 1rrnin iin iiiyyhf四四. 線性多步法線性多步法l多步法的概念多步法的概念單步法單步法: 計(jì)算計(jì)算yn+1時(shí)時(shí), 只用到前一步的信息只用到前一步的信息yn;RK方法也能提高精度方法也能提高精度, 但需要計(jì)算多個(gè)點(diǎn)的函數(shù)值但需要計(jì)算多個(gè)點(diǎn)的函數(shù)值, 計(jì)算量較大計(jì)算量較大.多步法多步法: 計(jì)算計(jì)算yn+1時(shí)

41、時(shí), 利用前幾步的已有信息利用前幾步的已有信息yn, yn-1, yn-2,從而從而提高計(jì)算精度提高計(jì)算精度多步法的基本思想多步法的基本思想.最常用的多步法最常用的多步法線性多步法的一般形式線性多步法的一般形式前幾步的信息前幾步的信息(已知已知)(,) (1, ,)kkkff xyknnnr其中其中,1,nnn rfff前幾步的前幾步的信息信息(已知已知)1nf當(dāng)前未知信息當(dāng)前未知信息如果如果 , 需用到需用到 fn+1 (未知未知), 則公式為則公式為隱式隱式的的;若若 , 則公式為則公式為顯式顯式的的.1010第第9章常微分方程數(shù)值解法章常微分方程數(shù)值解法復(fù)旦大學(xué)力學(xué)與工程科學(xué)系復(fù)旦大學(xué)力

42、學(xué)與工程科學(xué)系9-34u線性多步法的導(dǎo)出線性多步法的導(dǎo)出l基本方法基本方法將多步法公式在將多步法公式在xn處處Talor展開展開, 與與y(xn+1)在在xn處的處的Talor展開式相展開式相比比, 要求前面的若干項(xiàng)相同要求前面的若干項(xiàng)相同, 從而確定系數(shù)從而確定系數(shù)i ,i . .以線性兩步法為例以線性兩步法為例, r=1, 假設(shè)假設(shè)f(x)充分光滑充分光滑.( )( )()kknnyyx. (* *)l線性兩步法的推導(dǎo)過(guò)程線性兩步法的推導(dǎo)過(guò)程101111011()nnnnnnyyyhfff線性兩步法公式線性兩步法公式:Talor展開展開:()2( )()()()2!ppnnnnnnnyyy

43、 xyyxxxxxxp1()pnOxx其中其中,假設(shè)假設(shè), 都是精確的都是精確的. 由由(* *)式可得式可得,(),()(,) ()iiiiy xy xf x yin(4)(5)2345611()()23!4!5!nnnnnnnnyyyyyy xyy hhhhhO h第第9章常微分方程數(shù)值解法章常微分方程數(shù)值解法復(fù)旦大學(xué)力學(xué)與工程科學(xué)系復(fù)旦大學(xué)力學(xué)與工程科學(xué)系9-35l線性兩步法的推導(dǎo)過(guò)程線性兩步法的推導(dǎo)過(guò)程 (續(xù)續(xù)-1)1111(,)()nnnnff xyy x(4)(5)2345()2!3!4!nnnnnyyyyy hhhhO h(,)nnnnff xyy1111(,)()nnnnff

44、xyy x(4)(5)2345()2!3!4!nnnnnyyyyy hhhhO h101111011()nnnnnnyyyhfff代入兩步法公式代入兩步法公式:(4)(5)234510123!4!5!nnnnnnnnyyyyyyyy hhhhh(4)(5)234102!3!4!nnnnnnyyyhyy hhhhhy(4)(5)23412!3!4!nnnnnyyyhyy hhhh第第9章常微分方程數(shù)值解法章常微分方程數(shù)值解法復(fù)旦大學(xué)力學(xué)與工程科學(xué)系復(fù)旦大學(xué)力學(xué)與工程科學(xué)系9-36l線性兩步法的推導(dǎo)過(guò)程線性兩步法的推導(dǎo)過(guò)程 (續(xù)續(xù)-2)整理得整理得,y(xn+1)的的Talor展開式展開式:令令(

45、1)與與(2)式各項(xiàng)相同式各項(xiàng)相同, 得得:(4)(5)456()4!5!nnyyhhO h. (2)231()2!3!nnnnnyyy xyy hhh101()nnyy3111()622nyh(5)56111()()1202424nyhO h. (1)1101()nyh2111(2)nyh(4)41111()2466nyh選取選取 使其滿足使其滿足(3)式式前前p+1個(gè)個(gè)(p=1,2,3,4)方程方程,則可得則可得p階公式階公式.,ii . (3)01110111111111111122162261124662411111202424120第第9章常微分方程數(shù)值解法章常微分方程數(shù)值解法復(fù)旦大

46、學(xué)力學(xué)與工程科學(xué)系復(fù)旦大學(xué)力學(xué)與工程科學(xué)系9-37l線性兩步法的推導(dǎo)例子線性兩步法的推導(dǎo)例子(3)式滿足前式滿足前3式式.011101111111111111221622611246624二階公式（二階公式（p=2）01110111111122共共5個(gè)未知量個(gè)未知量, 有無(wú)窮個(gè)解有無(wú)窮個(gè)解. 取取011011,0,1 2,0滿足上面三個(gè)方程滿足上面三個(gè)方程. 二階公式二階公式:11()2nnnnhyyff梯形公式梯形公式隱式單步法隱式單步法四階公式（四階公式（p=4）(3)式滿足前式滿足前5式式.得唯一解得唯一解:1101 34 3,010,1,1111(4)3nnnnnhyyfff得四階公式

47、得四階公式,Simpson公式公式隱式二步法隱式二步法第第9章常微分方程數(shù)值解法章常微分方程數(shù)值解法復(fù)旦大學(xué)力學(xué)與工程科學(xué)系復(fù)旦大學(xué)力學(xué)與工程科學(xué)系9-38l一般形式線性多步法的系數(shù)確定一般形式線性多步法的系數(shù)確定令令(4)與與(5)式前式前p+1項(xiàng)相等項(xiàng)相等, 可得可得有有r+1個(gè)個(gè)i , r+2個(gè)個(gè)i , 共共2r+3個(gè)未知量個(gè)未知量. ()n inyy xih2(,)()()()2!nn in in innnyff xyy xihyyihih()(1)11()()()(1)!pppppnnyyihihO hpp()(1)212()()()()()2!(1)!pppppnnnnnyyyyy

48、ihihihihO hpp1=0= 1rrnin iin iiiyyhf. (4)()(1)121()()!(1)!pppppnnnnnyyy xyy hhhO hpp. (5)yn-i 及及fn-i 在在xn處作處作Talor展開展開, 得得第第9章常微分方程數(shù)值解法章常微分方程數(shù)值解法復(fù)旦大學(xué)力學(xué)與工程科學(xué)系復(fù)旦大學(xué)力學(xué)與工程科學(xué)系9-39l一般形式線性多步法的系數(shù)確定一般形式線性多步法的系數(shù)確定 (續(xù)續(xù))(1)(1)(1)11211()()()(1)!(1)!ppprrppppnnniiiiyyyhihhihO hppph0項(xiàng)項(xiàng)(即即yn項(xiàng)項(xiàng)) 的系數(shù)相等的系數(shù)相等=01riihk 項(xiàng)項(xiàng)

49、(即即 項(xiàng)項(xiàng)) 的系數(shù)相等的系數(shù)相等( )kny1=11()()1!(1)!kkrriiiiiikkk1=11()()1rrkkiiiiiki局部截?cái)嗾`差局部截?cái)嗾`差1nR1(1)121111()( +1)()()(1)!pprrpppnniiiihyRipiO hp共共2r+3個(gè)未知量個(gè)未知量,p+1個(gè)方程個(gè)方程.要求要求: p+12r+3由于由于p2r+2, 即即(4)式的線性多步法最多可達(dá)到式的線性多步法最多可達(dá)到2r+2階精度階精度.(1,2,)kp規(guī)定規(guī)定:001第第9章常微分方程數(shù)值解法章常微分方程數(shù)值解法復(fù)旦大學(xué)力學(xué)與工程科學(xué)系復(fù)旦大學(xué)力學(xué)與工程科學(xué)系9-40u常用的線性多步法常

50、用的線性多步法l阿達(dá)姆斯阿達(dá)姆斯 (Adams) 公式公式線性多步公式為線性多步公式為:1123(5559379)24nnnnnnhyyffff取取r=3, p=42r+2;12310,3=0331=111()()1iikkiiiiiki(1,2,)kp01即即,012312312312312(23)13(49)14(827)1310()1kiiki01235559,2424379,2424 四階四階Adams顯式公式顯式公式局部截?cái)嗾`差為局部截?cái)嗾`差為:5(5)335461111()5()()5!nniiiih yRiiO h5(5)6251()720nh yO h第第9章常微分方程數(shù)值解法

51、章常微分方程數(shù)值解法復(fù)旦大學(xué)力學(xué)與工程科學(xué)系復(fù)旦大學(xué)力學(xué)與工程科學(xué)系9-413=0331=111()()1iikkiiiiiki(1,2,)kpl阿達(dá)姆斯阿達(dá)姆斯 (Adams) 公式（續(xù)）公式（續(xù)）1112(9195)24nnnnnnhyyffff仍取仍取r=3, p=4; 并取并取12330,01即即,101211211211212(2)13(4)14(8)1211()1kiiki1012919,242451,2424 四階四階Adams隱式公式隱式公式注意注意:規(guī)定規(guī)定:001線性多步公式為線性多步公式為:局部截?cái)嗾`差為局部截?cái)嗾`差為:5(5)335461111()5()()5!nnii

52、iih yRiiO h5(5)619()720nh yO h 第第9章常微分方程數(shù)值解法章常微分方程數(shù)值解法復(fù)旦大學(xué)力學(xué)與工程科學(xué)系復(fù)旦大學(xué)力學(xué)與工程科學(xué)系9-423=0331=111()()1iikkiiiiiki(1,2,3,4)k l米爾恩米爾恩 (Milne) 公式公式13124(22)3nnnnnyyhfff取取r=3, p=4; 并取并取31即即,3130( 3)()1kkiiki012384,338,03 注意注意:規(guī)定規(guī)定:001線性多步公式為線性多步公式為:局部截?cái)嗾`差為局部截?cái)嗾`差為:5(5)335461111()5()()5!nniiiih yRiiO h5(5)614(

53、)45nh yO h01210,01231231231233192(23)1273(49)1814(827)1 Milne公式公式(四階顯式四階顯式)第第9章常微分方程數(shù)值解法章常微分方程數(shù)值解法復(fù)旦大學(xué)力學(xué)與工程科學(xué)系復(fù)旦大學(xué)力學(xué)與工程科學(xué)系9-432=0221=111()()1iikkiiiiiki(1,2,3,4)k l海明海明 (Hamming) 公式公式121113(9)(2)88nnnnnnyyyh fff取取r=2, p=40, 使得使得:整體截?cái)嗾`差為整體截?cái)嗾`差為,1111()()(, )nnnnnnney xyy xyhxyh1()()(, (), )nnnny xy xh

54、xy xh () (, (), )(, )nnnnnny xyhxy xhxyh11()(, (), )(, )nnnnnnnneRy xyhxy xhxyh1nR設(shè)顯式單步法具有設(shè)顯式單步法具有p階精度階精度, 其增量函數(shù)其增量函數(shù) 關(guān)于關(guān)于y滿足滿足Lipschitz 條件條件, 初值初值y0=y(x0)是精確的是精確的. 則顯式單步法的整體截?cái)嗾`差為則顯式單步法的整體截?cái)嗾`差為:ne1212( , )( , )x y hx y hL yy利用利用Lipschitz 條件條件:L 為為常數(shù)常數(shù).111() ()(, (), )pnnnnnRy xy xhxy xhMh (2)(0,1,1)

55、nN，第第9章常微分方程數(shù)值解法章常微分方程數(shù)值解法復(fù)旦大學(xué)力學(xué)與工程科學(xué)系復(fù)旦大學(xué)力學(xué)與工程科學(xué)系9-61定理定理9.1 證明證明 (續(xù)續(xù))11()nnnnneRehL y xy1(1)nnRehL同理同理,11(1),pnneMhehL1111(1)1(1)npnhLeMhhL11(1)pnneMhehL1nn+1101(1)(1) (1)pneMhhLhLehL00()0y xy因?yàn)樵谝驗(yàn)樵赼,b中求解中求解, 故有故有:(1)()nhba(1)banh利用不等式利用不等式:1,(0)xxex()b ahLL b ahee1(1)(1)b anhhLhL(取取x=hL)1(1)1nphL

56、MhL()11()L b appnMeehO hL證畢證畢 #第第9章常微分方程數(shù)值解法章常微分方程數(shù)值解法復(fù)旦大學(xué)力學(xué)與工程科學(xué)系復(fù)旦大學(xué)力學(xué)與工程科學(xué)系9-62推論推論1連續(xù)連續(xù), 且關(guān)于且關(guān)于y滿足滿足Lipschitz 條件條件, 則顯式單步法是收斂的則顯式單步法是收斂的.設(shè)顯式單步法具有設(shè)顯式單步法具有p階精度階精度, 其增量函數(shù)其增量函數(shù) 在區(qū)域在區(qū)域( , , )x y h0:, 0Saxbyhh 證明證明: 由定理由定理9.1的結(jié)論可直接推得的結(jié)論可直接推得.由此可得由此可得, 當(dāng)當(dāng)f(x,y)在區(qū)域在區(qū)域D: axb, -y+上連續(xù)上連續(xù), 且關(guān)且關(guān)于于y滿足滿足Lipsch

57、itz 條件時(shí)條件時(shí):Euler法、改進(jìn)法、改進(jìn)Euler法及各階法及各階RK方法的增量函數(shù)方法的增量函數(shù) 在在區(qū)域區(qū)域S上連續(xù)上連續(xù), 且關(guān)于且關(guān)于y滿足滿足Lipschitz 條件條件, 因此都是收斂的因此都是收斂的.( , , )x y h1=1pnniiiyyhc K例如例如, RK方法方法, 11=1(,),(,)innininijjjKf xyKf xa h yhb K其中其中,=1( , , )piiix y hc K由由1212( ,)( ,)f x yf x yL yy1212( , )( , )( )x y hx yhg Lyy其中其中g(shù)(L)是與是與L有關(guān)的常數(shù)有關(guān)的常數(shù)

58、. 故滿足故滿足Lipschitz條件條件.第第9章常微分方程數(shù)值解法章常微分方程數(shù)值解法復(fù)旦大學(xué)力學(xué)與工程科學(xué)系復(fù)旦大學(xué)力學(xué)與工程科學(xué)系9-63l收斂性的判斷（續(xù)收斂性的判斷（續(xù)-2）定理定理9.2證明證明: (略略)設(shè)增量函數(shù)設(shè)增量函數(shù) 在區(qū)域在區(qū)域S上連續(xù)上連續(xù), 且關(guān)于且關(guān)于y滿足滿足Lipschitz 條條件件, 則顯式單步法收斂的充分必要條件是則顯式單步法收斂的充分必要條件是, 滿足相容性條件滿足相容性條件:( , , )x y h( , ,0)( , )x yf x y第第9章常微分方程數(shù)值解法章常微分方程數(shù)值解法復(fù)旦大學(xué)力學(xué)與工程科學(xué)系復(fù)旦大學(xué)力學(xué)與工程科學(xué)系9-64u穩(wěn)定性穩(wěn)

59、定性l穩(wěn)定性的概念穩(wěn)定性的概念收斂性收斂性: 只考慮數(shù)值方法的整體截?cái)嗾`差只考慮數(shù)值方法的整體截?cái)嗾`差, 不考慮計(jì)算過(guò)程中的不考慮計(jì)算過(guò)程中的舍入誤差舍入誤差. 即在假設(shè)每一步計(jì)算結(jié)果都是精確的前提下即在假設(shè)每一步計(jì)算結(jié)果都是精確的前提下, 討論當(dāng)討論當(dāng)h0時(shí)時(shí), 數(shù)值解是否趨于精確數(shù)值解是否趨于精確解解.特別地特別地, 如果絕對(duì)穩(wěn)定區(qū)域包含整個(gè)左半平面如果絕對(duì)穩(wěn)定區(qū)域包含整個(gè)左半平面, 則稱該方法是則稱該方法是A-穩(wěn)穩(wěn)定的定的, 絕對(duì)穩(wěn)定區(qū)域與實(shí)軸的交稱為絕對(duì)穩(wěn)定區(qū)間絕對(duì)穩(wěn)定區(qū)域與實(shí)軸的交稱為絕對(duì)穩(wěn)定區(qū)間.穩(wěn)定性穩(wěn)定性: 討論當(dāng)某一節(jié)點(diǎn)上的計(jì)算值討論當(dāng)某一節(jié)點(diǎn)上的計(jì)算值yn存在一個(gè)誤差存在一

60、個(gè)誤差(擾動(dòng)擾動(dòng)), 則則該擾動(dòng)是否會(huì)隨著該擾動(dòng)是否會(huì)隨著n的增加而不斷增加的增加而不斷增加.l穩(wěn)定性的定義穩(wěn)定性的定義定義定義: 對(duì)于給定的微分方程和給定的步長(zhǎng)對(duì)于給定的微分方程和給定的步長(zhǎng)h, 某種方法計(jì)算節(jié)點(diǎn)某種方法計(jì)算節(jié)點(diǎn)值值yn時(shí)產(chǎn)生誤差為時(shí)產(chǎn)生誤差為, 由此所引起后續(xù)節(jié)點(diǎn)值由此所引起后續(xù)節(jié)點(diǎn)值ym的誤差的誤差m, 若若 , 則稱該方法對(duì)所用步長(zhǎng)則稱該方法對(duì)所用步長(zhǎng)h是是絕對(duì)穩(wěn)定絕對(duì)穩(wěn)定的的. 使得方法絕對(duì)使得方法絕對(duì)穩(wěn)定的穩(wěn)定的h及方法中的參數(shù)全體稱為該方法的及方法中的參數(shù)全體稱為該方法的絕對(duì)穩(wěn)定區(qū)域絕對(duì)穩(wěn)定區(qū)域.m第第9章常微分方程數(shù)值解法章常微分方程數(shù)值解法復(fù)旦大學(xué)力學(xué)與工程科