數值計算方法-第9章-常微分方程數值解

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

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

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

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

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

6、常微分方程數值解法章常微分方程數值解法復旦大學力學與工程科學系復旦大學力學與工程科學系9-5l用數值積分方法用數值積分方法用左端點的矩形公式計算積分用左端點的矩形公式計算積分11()()( , )nnxnnxy xy xf x y dx( , )dyf x ydx兩端兩端積分積分nx1nx( )f x1(,)nnnnyyhf xy與差商近似的公式相同與差商近似的公式相同得得y(xn+1)的近似值的近似值yn+1用梯形公式計算積分用梯形公式計算積分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章常微分方程數值解法章常微分方程數值解法復旦大學力學與工程科學系復旦大學力學與工程科學系9-6l用用Talor 展開近似的方法展開近似的方法用用y(xn)的近似值的近似值yn代入上式代入上式, 可得可得y(xn+1)的近似值的近似值yn+1:2()11()( )( )( )( )2!ppy xhy xhy xh yxh yxp1(, )nnnnyyhxyh上面三種方法都能導出常微分方程數值解的公式上面三種方法都能導出常微分方程數值解的公式, 其中其中

8、Talor 展展開方法還可以用來估計截斷誤差開方法還可以用來估計截斷誤差, 所以推導時都用該方法所以推導時都用該方法.與差商近似的公式相同與差商近似的公式相同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局部截斷誤差或局部截斷誤差或簡稱截斷誤差簡稱截斷誤差第第9章常微分方程數值解法章常微分

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

10、的 幾 何 解 釋P0P1P2P3P4Pn( )yf xOx0 x1 x2 x3 x4 xnxyEuler 方法是顯式的方法是顯式的, 可直接遞可直接遞推求解推求解.第第9章常微分方程數值解法章常微分方程數值解法復旦大學力學與工程科學系復旦大學力學與工程科學系9-8l向后向后Euler 方法方法向后向后Euler方法方法這是隱式公式這是隱式公式, 一般用迭代法求解一般用迭代法求解:若用向后差商代替導數若用向后差商代替導數, 即即: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 方法方法若用中心差商代替導數若用中心差商代替導數, 即即:211()()()()2nnny xy xy xO hh中點中點Euler方法方法3112()()nnnyyhy xO h112(,)nnnnyyhf xy因為要用到前面兩步的結果因為要用到前面兩步的結果yn-1,yn, 故又稱為故又稱為Euler兩步公式兩步公式.第第9章常微分方程數值解法章常微分方程數值解法復旦大學力學與工程科學系復旦大學力學與工程

12、科學系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, 使得使得整體截斷誤差整體截斷誤差: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ù)函數由閉區(qū)間上連續(xù)函數的性質

13、的性質, Rn+1在在a,b區(qū)間有界區(qū)間有界.第第9章常微分方程數值解法章常微分方程數值解法復旦大學力學與工程科學系復旦大學力學與工程科學系9-14lEuler 方法的整體截斷誤差方法的整體截斷誤差 (續(xù)續(xù))反復遞推得反復遞推得,211(1)2nneh MhL e221111(1)(1)22nneh MhLh MhL e21001(1)(1)2nknkh MhLhLe021(1)12(1)1nh MhLhL1(1)12nhMhLL因為在因為在a,b中求解中求解, 故有故有:(1)()nhba(1)banh1(1)(1)b anhhLhL利用不等式利用不等式:1,(0)xxex()b ahLL

14、b ahee(取取x=hL)()112L b anhMeeLEuler方法的整體截斷誤差與方法的整體截斷誤差與h同階同階. 當當h0時時, eN0第第9章常微分方程數值解法章常微分方程數值解法復旦大學力學與工程科學系復旦大學力學與工程科學系9-15二二. 改進的歐拉改進的歐拉 (Euler) 方法方法u 梯形公式梯形公式l梯形公式的導出梯形公式的導出梯形公式梯形公式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章常微分方程數值解法章常微分方程數值解法復旦大學力學與工程科學系復旦大學力學與工程科學系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)關于關于y滿足滿足Lipschitz條件條件, (1)( )( )(1)111111(,)(,)2kkkknnnnnnhyyf xyf xy( )(1)112kknnhL yy(1)(0)

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

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

19、3.2置置xx+h; 計算計算yq=y+hf(x,yp);3.3置置y0.5(yp+yq);3.4輸出輸出(x,y);End for il改進改進Euler法的誤差法的誤差可以證明可以證明, 改進改進Euler 法的局部截斷誤差為法的局部截斷誤差為O(h3), 故屬于二階方故屬于二階方法法.第第9章常微分方程數值解法章常微分方程數值解法復旦大學力學與工程科學系復旦大學力學與工程科學系9-19l改進改進Euler 法的算例法的算例從結果可以看出從結果可以看出: 精確度比精確度比Euler法明顯提高法明顯提高.誤誤差差越越來來越越大大取步長取步長h=0.1; 結果如下結果如下: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章常微分方程數值解法章常微分方程數值解法復旦大學力學與工程科學系復旦大學力學與工程科學系9-20三三. 龍格龍格-庫塔庫塔 (Runge-Kutta) 法法u Runge-Kutta 法的基本思想法的基本思想l如何提高精度如何提高精度局部截斷誤差由局部截斷誤差由Talor展開式的余項決定展開式的余項決定, 一般地一般地,2()1()()()2!ppnnnnnhhyyhy xyxyxp具有具有p 階精度階精度, 局部截斷誤差為局部截斷誤差為:1()pO h其中其中, ( )( , )( , ),xyyxfx yfx y y( )( , ),

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

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

24、 )f x y第第9章常微分方程數值解法章常微分方程數值解法復旦大學力學與工程科學系復旦大學力學與工程科學系9-22uRKRK方法的構造方法的構造l構造方法構造方法1=111=1(,)(,)pnniiinniininijjjyyhc KKf xyKf xa h yhb K其其中中，一般地一般地, 設設RK近似公式為近似公式為:其中其中, 為待定參數為待定參數.,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章常微分方程數值解法章常微分方程數值解法復旦大學力學與工程科學系復旦大學力學與工程科學系9-23l二階二階RK方法（續(xù)）方法（續(xù)）理論上可以證明理論上可以證明, 無論怎樣選取參數無論怎樣選取參數c1, c2, a

26、2, b21, 上面的公式不可上面的公式不可能有更高的精度能有更高的精度. 即即, 每步計算兩個函數值每步計算兩個函數值, 只能得出二階公式只能得出二階公式.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個未知量個未知量, 有無窮多個解有無窮多個解.代入代入RK近似公式近似公式, 其局部截斷誤差其局部截斷誤差都是都是O(h3). 統(tǒng)稱為統(tǒng)稱為二階方法二階方法.21222121 (2)1bacacc 第第9章常微分方程數值解法章常微分方程數值解法復旦大學力學與工程科學系復旦大學力學與工程科學系9-24l幾個常用的二階幾個常用的二階RK公式公式122222111 21 2ccc ac b改進改進Euler公式公式取取122211 2,1ccab11122122211()(,)(,)nnnnnnyyhc Kc KKf xyKf

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

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

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

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

32、科學系復旦大學力學與工程科學系9-28l經典四階經典四階Kutta 方法的算法描述方法的算法描述1.輸入輸入a,b, f(x,y), 區(qū)間等分數區(qū)間等分數N, 初值初值y0;3.for i=1 to N step 1 計算各點的函數值計算各點的函數值4.停機停機.2.計算步長計算步長 h=(b-a)/N; 置置xa, yy0; 輸出輸出(x,y);3.1 計算計算 K1= f(x,y);3.2置置xx+0.5*h; 計算計算:K2= f(x,y+0.5*h*K1);K3= f(x,y+0.5*h*K2); 3.3置置x a+i*h; 計算計算: 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 計算各點的函數值計算各點的函數值第第9章常微分方程數值解法章常微分方程數值解法復旦大學力學與工程科學系復旦大學力學與工程科學系9-29l經典四階經典四階RK RK 公式算例公式算例從結果可以看出從結果可以看出: 精確度比精確度比改進改進Euler法明顯提高法明顯提高.誤誤差差越越來來越越大大取步長取步長h=0.1; 結果如下結果如下: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章常微分方程數值解法章常微

35、分方程數值解法復旦大學力學與工程科學系復旦大學力學與工程科學系9-30u變步長的變步長的RKRK方法方法l為何要變步長為何要變步長由于由于y(x)是不均勻的是不均勻的, 采用等不長計算時采用等不長計算時, 有些點上精度較高有些點上精度較高, 而而有些點則精度較低有些點則精度較低. 是否能根據精度要求是否能根據精度要求, 自動調整當前計算步的自動調整當前計算步的步長呢步長呢?局部截斷誤差局部截斷誤差O(hp+1), 即即( )111()hpnny xychnyp階方法階方法,步長步長h( )1hnyc與與 有關有關, 假設假設 內變內變化不大化不大, 則可看作常數則可看作常數.(1)1( ),p

36、nnyxx(1)1( ),pnnyxx在在. (1)nyp階方法階方法,步長步長h/2(2)1 2hnyp階方法階方法,步長步長h/2(2)1hny截斷誤差截斷誤差O(h/2)p+1, 即即1( /2)11()22phnnhy xyc. (2)l從誤差的分析說起從誤差的分析說起如果分兩步推進如果分兩步推進,第第9章常微分方程數值解法章常微分方程數值解法復旦大學力學與工程科學系復旦大學力學與工程科學系9-31l變步長控制方法變步長控制方法-1-1具體實施步驟具體實施步驟( )( /2)( )1111(21) ()22phphphnnnny xyyy( )( /2)( )11112()21phhh

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

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

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

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

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

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

43、 xyyxxxxxxp1()pnOxx其中其中,假設假設, 都是精確的都是精確的. 由由(* *)式可得式可得,(),()(,) ()iiiiy xy xf x yin(4)(5)2345611()()23!4!5!nnnnnnnnyyyyyy xyy hhhhhO h第第9章常微分方程數值解法章常微分方程數值解法復旦大學力學與工程科學系復旦大學力學與工程科學系9-35l線性兩步法的推導過程線性兩步法的推導過程 (續(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章常微分方程數值解法章常微分方程數值解法復旦大學力學與工程科學系復旦大學力學與工程科學系9-36l線性兩步法的推導過程線性兩步法的推導過程 (續(xù)續(xù)-2)整理得整理得,y(xn+1)的的Talor展開式展開式:令令(

45、1)與與(2)式各項相同式各項相同, 得得:(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個個(p=1,2,3,4)方程方程,則可得則可得p階公式階公式.,ii . (3)01110111111111111122162261124662411111202424120第第9章常微分方程數值解法章常微分方程數值解法復旦大

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

47、得四階公式,Simpson公式公式隱式二步法隱式二步法第第9章常微分方程數值解法章常微分方程數值解法復旦大學力學與工程科學系復旦大學力學與工程科學系9-38l一般形式線性多步法的系數確定一般形式線性多步法的系數確定令令(4)與與(5)式前式前p+1項相等項相等, 可得可得有有r+1個個i , r+2個個i , 共共2r+3個未知量個未知量. ()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章常微分方程數值解法章常微分方程數值解法復旦大學力學與工程科學系復旦大學力學與工程科學系9-39l一般形式線性多步法的系數確定一般形式線性多步法的系數確定 (續(xù)續(xù))(1)(1)(1)11211()()()(1)!(1)!ppprrppppnnniiiiyyyhihhihO hppph0項項(即即yn項項) 的系數相等的系數相等=01riihk 項項

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

50、用的線性多步法l阿達姆斯阿達姆斯 (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顯式公式顯式公式局部截斷誤差為局部截斷誤差為:5(5)335461111()5()()5!nniiiih yRiiO h5(5)6251()720nh yO h第第9章常微分方程數值解法

51、章常微分方程數值解法復旦大學力學與工程科學系復旦大學力學與工程科學系9-413=0331=111()()1iikkiiiiiki(1,2,)kpl阿達姆斯阿達姆斯 (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線性多步公式為線性多步公式為:局部截斷誤差為局部截斷誤差為:5(5)335461111()5()()5!nnii

52、iih yRiiO h5(5)619()720nh yO h 第第9章常微分方程數值解法章常微分方程數值解法復旦大學力學與工程科學系復旦大學力學與工程科學系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線性多步公式為線性多步公式為:局部截斷誤差為局部截斷誤差為:5(5)335461111()5()()5!nniiiih yRiiO h5(5)614(

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

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

55、nN，第第9章常微分方程數值解法章常微分方程數值解法復旦大學力學與工程科學系復旦大學力學與工程科學系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因為在因為在a,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章常微分方程數值解法章常微分方程數值解法復旦大學力學與工程科學系復旦大學力學與工程科學系9-62推論推論1連續(xù)連續(xù), 且關于且關于y滿足滿足Lipschitz 條件條件, 則顯式單步法是收斂的則顯式單步法是收斂的.設顯式單步法具有設顯式單步法具有p階精度階精度, 其增量函數其增量函數 在區(qū)域在區(qū)域( , , )x y h0:, 0Saxbyhh 證明證明: 由定理由定理9.1的結論可直接推得的結論可直接推得.由此可得由此可得, 當當f(x,y)在區(qū)域在區(qū)域D: axb, -y+上連續(xù)上連續(xù), 且關且關于于y滿足滿足Lipsch

57、itz 條件時條件時:Euler法、改進法、改進Euler法及各階法及各階RK方法的增量函數方法的增量函數 在在區(qū)域區(qū)域S上連續(xù)上連續(xù), 且關于且關于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(L)是與是與L有關的常數有關的常數

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

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

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