《計算方法第八章》ppt課件

1、第八章第八章 常微分方程初值常微分方程初值問題數(shù)值解法問題數(shù)值解法本章主要研討常微分方程初值問題的數(shù)值求解：本章主要研討常微分方程初值問題的數(shù)值求解：通常，假設函數(shù)通常，假設函數(shù) f 關(guān)于第二個變量滿足李普希茨條件關(guān)于第二個變量滿足李普希茨條件L條條件，即為存在常數(shù)件，即為存在常數(shù) L 0，使得，使得0( )( , ( )( )y xf x y xaxby ay1212( ,)( ,)f x yf x yL yy第一節(jié)第一節(jié) 普通概念普通概念1.1 歐拉法及其簡單改良歐拉法及其簡單改良0( )( , ( )( )y xf x y xaxby ay方法：選擇適當?shù)墓?jié)點，用差分近似微分，將方程離散

2、化，方法：選擇適當?shù)墓?jié)點，用差分近似微分，將方程離散化，從而求在這些節(jié)點上的函數(shù)值。從而求在這些節(jié)點上的函數(shù)值。012111NNnnnnnaxxxxxbhxxxx稱為到的步長（通常取為常數(shù)h）1()()(, ()(,)nnnnnny xy xf xy xf xyh歐拉歐拉方法方法100(,),()nnnnyyhf xyy xy例子：例子：(01)(0)1yyxy 10( , )(1)1nnnnf x yyyyhyh yy xye精確解為下面我們分別取步長為下面我們分別取步長為0.1與與0.01進展計算，進展計算，計算結(jié)果顯示在下面的圖中。計算結(jié)果顯示在下面的圖中。步長為步長為0.1的計算結(jié)果。

3、的計算結(jié)果。步長為步長為0.01的計算結(jié)果的計算結(jié)果0.01 0.99005 0.99 0.1 0.90484 0.90438 0.2 0.81873 0.81791 0.3 0.74082 0.7397 0.4 0.67032 0.66897 0.41 0.66365 0.66228 0.59 0.55433 0.55268 0.6 0.54881 0.54716 0.9 0.40657 0.40473 0.91 0.40252 0.40068 0.99 0.37158 0.36973 1 0.36788 0.36603 DOUBLE PRECISION h,y(0：100) OPEN(20

4、,FILE=OUTPUT.DAT,STATUS=UNKNOWN)h=1.0/100y(0)=1.0do 10 i=1,100 y(i)=y(i-1)*(1.0-h)write(20,*) i*h,y(i)10 continueEND精度更高的歐拉公式：精度更高的歐拉公式：方法：選擇計算導數(shù)值方法：選擇計算導數(shù)值的精度更好的差分格式。的精度更好的差分格式。11()2()()2nnnnnyyy xhy xhy xhh3222()( )/3!( )()()()212nnny x hyhyy xhy xO hh11002(,),()nnnnyyhf xyy xy歐拉中點公式歐拉中點公式利用中點公式求解

5、微分方程時，有一個問題，就是計算時需利用中點公式求解微分方程時，有一個問題，就是計算時需求兩個迭代初值！求兩個迭代初值！對于這個問題，我們可以先用歐拉公式，經(jīng)過給定的初值計對于這個問題，我們可以先用歐拉公式，經(jīng)過給定的初值計算出第一個點的值，然后在利用這兩個第一和第二個點算出第一個點的值，然后在利用這兩個第一和第二個點的值進展計算，直到計算出全部節(jié)點上的值。的值進展計算，直到計算出全部節(jié)點上的值。下面，我們用中點公式求解剛剛的例子，計算的步長取下面，我們用中點公式求解剛剛的例子，計算的步長取0.01，可以看到，計算的精度比較高?？梢钥吹?，計算的精度比較高。此時，計算公式為此時，計算公式為11(

6、 , )2nnnf x yyyyhy 100(1),1yh yy計計算算結(jié)結(jié)果果0.02 0.9802 0.98020.1 0.90484 0.904840.2 0.81873 0.818740.3 0.74082 0.740840.4 0.67032 0.670350.5 0.60653 0.606560.6 0.54881 0.548850.7 0.49659 0.496630.8 0.44933 0.449380.9 0.40657 0.406631.2 歐拉方法的其他改良歐拉方法的其他改良微分方程數(shù)值解的關(guān)鍵在于對導數(shù)的處置，可以用差分來近似微分方程數(shù)值解的關(guān)鍵在于對導數(shù)的處置，可以用

7、差分來近似導數(shù)，也可以經(jīng)過積分，將導數(shù)項化掉。導數(shù)，也可以經(jīng)過積分，將導數(shù)項化掉。對于方程：對于方程：首先，作出劃分首先，作出劃分設曾經(jīng)求出第設曾經(jīng)求出第 n 個節(jié)點的函數(shù)值個節(jié)點的函數(shù)值 ，在區(qū)間，在區(qū)間 上對上對方程兩邊積分方程兩邊積分容易看出，要求第容易看出，要求第 n+1 個節(jié)點的函數(shù)值，關(guān)鍵在于選擇適當個節(jié)點的函數(shù)值，關(guān)鍵在于選擇適當?shù)姆e分公式計算積分！的積分公式計算積分！( )( , ( )y xf x y x0121NNaxxxxxb1,nnxxny11()()( , ( )nnxnnxy xy xf x y x dx1如選擇下矩形公式，那么如選擇下矩形公式，那么這正是前面的歐拉

8、公式。這正是前面的歐拉公式。1(,)nnnnyyf xy h111(,)nnnnyyf xyh111 (,)(,)/2nnnnnnyyf xyf xyh2如選擇上矩形公式，那么如選擇上矩形公式，那么這是所謂的后退歐拉公式。這是所謂的后退歐拉公式。3如選擇梯形公式，那么如選擇梯形公式，那么這是所謂的歐拉梯形公式。這是所謂的歐拉梯形公式。直接利用曾經(jīng)求得的知節(jié)點上的值計算未知節(jié)點上的函數(shù)值直接利用曾經(jīng)求得的知節(jié)點上的值計算未知節(jié)點上的函數(shù)值的算法稱為顯式法。的算法稱為顯式法。例如：歐拉公式、歐拉中點公式例如：歐拉公式、歐拉中點公式計算未知節(jié)點上的函數(shù)值時，用到了未知節(jié)點上的函數(shù)值，計算未知節(jié)點上的

9、函數(shù)值時，用到了未知節(jié)點上的函數(shù)值，這種算法稱為隱式法。這種算法稱為隱式法。例如：后退歐拉法、歐拉梯形公式例如：后退歐拉法、歐拉梯形公式顯然，利用隱式法求微分方程的數(shù)值解是，需求從表達式中顯然，利用隱式法求微分方程的數(shù)值解是，需求從表達式中反解未知節(jié)點上的函數(shù)值。反解未知節(jié)點上的函數(shù)值。1.3 隱式法的詳細計算：隱式法的詳細計算：例如歐拉梯形公式例如歐拉梯形公式用迭代法計算用迭代法計算n+1步的值。步的值。1簡單迭代簡單迭代收斂的條件：收斂的條件：11111 (,)(,)/2(,)/2(,)/2nnnnnnnnnnnyyf xyf xyhyf xyhf xyh(0)1(,)nnnnyyhf x

10、y(1)( )111(,)/2(,)/2kknnnnnnyyf xyhf xyh/21Lh2牛頓迭代牛頓迭代必需指出，在真正計算中，常用的是簡單迭代，而且只迭代一必需指出，在真正計算中，常用的是簡單迭代，而且只迭代一步，迭代初值步，迭代初值 稱為預測，迭代稱為校正，這樣組成的稱為預測，迭代稱為校正，這樣組成的一組計算公式稱為預測校正公式。一組計算公式稱為預測校正公式。( )( )(1)( )11111( )11(,)/2(,)/21(,)/2kkkknnnnnnnnkynnyyf xyhf xyhyyfxyh(0)1ny(0)1(0)111(,) (,)(,)/2nnnnnnnnnnyyhf

11、xyyyf xyf xyh預測校正公式也稱為改良的歐拉法，將上面的組合公式改預測校正公式也稱為改良的歐拉法，將上面的組合公式改寫為：寫為：留意到留意到 ，將上式進一步改寫為：，將上式進一步改寫為：這是我們最終運用的計算格式。這是我們最終運用的計算格式。11 (,)(,(,)/2nnnnnnnnyyf xyf xyhf xyh1121211()2(,)(,)nnnnnnyyKKKhf xyKhf xh yK1nnxxh例子：例子：取步長為取步長為0.1計算，結(jié)果如圖。計算，結(jié)果如圖。(01)(0)1yyxy 11212111()2(,)(,)()nnnnnnnnyyKKKhf xyhyKhf x

12、h yKh yK 圖：圖： DOUBLE PRECISION h,y(0:10),ak1,ak2 OPEN(20,FILE=OUTPUT1.DAT,STATUS=UNKNOWN)h=1.0/10y(0)=1.0do 10 i=1,10 ak1=-h*y(i-1) ak2=-h*(y(i-1)+ak1) y(i)=y(i-1)+(ak1+ak2)/2.010 continue do 20 i=0,10 write(20,*) i*h,y(i),exp(-i*h)20 continueEND同理，對于后退歐拉公式同理，對于后退歐拉公式有預測校正公式有預測校正公式或改寫為：或改寫為：111(,)nn

13、nnyyf xyh(0)1(0)111(,)(,)nnnnnnnnyyhf xyyyf xyh11(,)(,)nnnnnnKhf xyyyf xyK h用此法解前面的例子用此法解前面的例子步長步長0.1步長步長0.011.4 誤差估計誤差估計定義：利用第定義：利用第n個節(jié)點的函數(shù)值，經(jīng)過數(shù)值方法計算第個節(jié)點的函數(shù)值，經(jīng)過數(shù)值方法計算第n+1個個節(jié)點的函數(shù)近似值，所引起的誤差，稱為節(jié)點的函數(shù)近似值，所引起的誤差，稱為n+1個節(jié)點上的部個節(jié)點上的部分截斷誤差。分截斷誤差。我們記我們記 為為n+1個節(jié)點上函數(shù)的準確值，個節(jié)點上函數(shù)的準確值， 為數(shù)值為數(shù)值解，解，那么部分截斷誤差為：那么部分截斷誤差為

14、：如部分截斷誤差為如部分截斷誤差為 ，稱為具有，稱為具有p階部分截斷誤差。階部分截斷誤差。1ny1()ny x111()nnny xy()pO h歐拉方法的誤差分析：歐拉方法的誤差分析：12()()()()()()( )/2()()( ) /2nnnnnnnny xy xy xhy xhhy xy x hyhy xy xyhh1()()()( )nnny xy xy xO hh121()()()( )(, ()()()(, ()()nnnnnnnnny xy xy xO hf xy xhy xy xhf xy xO h省略余項得公式：省略余項得公式：1()(, ()nnnnyy xhf xy

15、x2111()()nnny xyO h完全類似的可以得到完全類似的可以得到后退歐拉公式的部分截斷誤差為：后退歐拉公式的部分截斷誤差為：歐拉中點公式的部分截斷誤差為：歐拉中點公式的部分截斷誤差為：歐拉梯形公式的部分截斷誤差為：歐拉梯形公式的部分截斷誤差為：2111()()nnny xyO h3111()()nnny xyO h3111()()nnny xyO h定義：由初值引起的第定義：由初值引起的第 n 個節(jié)點的近似值與準確值之間的誤差個節(jié)點的近似值與準確值之間的誤差稱為第稱為第 n 個節(jié)點整體誤差。個節(jié)點整體誤差。定理：設下面求解微分方程的數(shù)值計算方法部分截斷誤差為定理：設下面求解微分方程的

16、數(shù)值計算方法部分截斷誤差為p+1階，且函數(shù)階，且函數(shù) 關(guān)于關(guān)于 y 滿足利普希茨條件，滿足利普希茨條件，同時初值是準確的，那么整體截斷誤差為同時初值是準確的，那么整體截斷誤差為p階。階。歐拉公式、后退歐拉公式的整體截斷誤差為歐拉公式、后退歐拉公式的整體截斷誤差為 1 階。階。歐拉中點公式、歐拉梯形公式的整體截斷誤差為歐拉中點公式、歐拉梯形公式的整體截斷誤差為 2 階。階。1(, )nnnnyyhxy h( , , )x y h( , , )( , , )x y hx y hL yy微分方程數(shù)值解法的進一步改良。再回到恒等式微分方程數(shù)值解法的進一步改良。再回到恒等式假設取假設取 作為節(jié)點，將被積

17、函數(shù)用插值多作為節(jié)點，將被積函數(shù)用插值多項式來近似，用插值多項式帶到積分中去求出積分，那么可以項式來近似，用插值多項式帶到積分中去求出積分，那么可以得得到所謂的亞當斯到所謂的亞當斯(Adams)顯式公式顯式公式部分截斷誤差：部分截斷誤差：11()( )( , ( )iixiixy xy xf x y x dx1,iii kx xx1011()iiiiki khyyb fb fb fA1(1) ( )kkkiR yB hy類似地，假設取類似地，假設取 作為作為節(jié)點，可得亞當斯節(jié)點，可得亞當斯(Adams)隱式公式隱式公式部分截斷誤差：部分截斷誤差：1,i iiii kxx xx*10112()i

18、iiiiki khyyb fb fb fb fA*2(2)1 ( )kkkiR yBhy進一步，假設我們將恒等式中的積分區(qū)間改為進一步，假設我們將恒等式中的積分區(qū)間改為 ，并在此區(qū)間上用辛普森公式，可得，并在此區(qū)間上用辛普森公式，可得1,i iixx1111()()( , ( )iixiixy xy xf x y x dx111143iiiiihyyfff1.5 絕對穩(wěn)定性絕對穩(wěn)定性一個常微分方程數(shù)值解法稱為絕對穩(wěn)定，是指在某一步如第一個常微分方程數(shù)值解法稱為絕對穩(wěn)定，是指在某一步如第一步產(chǎn)生的誤差如計算機的存儲誤差，在計算中會逐一步產(chǎn)生的誤差如計算機的存儲誤差，在計算中會逐漸減小。漸減小。稱

19、方程稱方程 為實驗方程，設計算步長為為實驗方程，設計算步長為 h ，設在計算開，設在計算開場時產(chǎn)生誤差存儲誤差，此誤差在以后會逐漸減弱，我場時產(chǎn)生誤差存儲誤差，此誤差在以后會逐漸減弱，我們稱該算法相對于們稱該算法相對于 是絕對穩(wěn)定的，是絕對穩(wěn)定的， 的全體的全體稱為該算法的絕對穩(wěn)定域。稱為該算法的絕對穩(wěn)定域。 yy hhh歐拉法的絕對穩(wěn)定區(qū)域歐拉法的絕對穩(wěn)定區(qū)域 后退歐拉法的絕對穩(wěn)定區(qū)域后退歐拉法的絕對穩(wěn)定區(qū)域11h11h1.6 部分截斷誤差的適用估計部分截斷誤差的適用估計1用兩種階數(shù)一樣的算法求解，計算出用兩種階數(shù)一樣的算法求解，計算出n+1步的近似值，步的近似值，從而得到部分截斷誤差估計。

20、從而得到部分截斷誤差估計。2用同樣的公式，用不同步長計算出用同樣的公式，用不同步長計算出n+1步的近似值，從步的近似值，從而得到部分截斷誤差估計。而得到部分截斷誤差估計。1.7 隱式法隱式法隱式法具有較好的絕對穩(wěn)定性！隱式法具有較好的絕對穩(wěn)定性！只不過在運用隱式法的時候，需求進展迭代，或者運用預測只不過在運用隱式法的時候，需求進展迭代，或者運用預測校正計算格式。校正計算格式。第二節(jié)第二節(jié) 泰勒級數(shù)法與龍格庫塔法泰勒級數(shù)法與龍格庫塔法對于方程：對于方程：取計算步長為取計算步長為 h ，那么，那么 ，將函數(shù)進展泰勒展，將函數(shù)進展泰勒展開開如函數(shù)如函數(shù) y( x )有有p+1階導數(shù)，容易得到階導數(shù)，

21、容易得到p階泰勒級數(shù)展開法：階泰勒級數(shù)展開法：0( )( , ( )( )y xf x y xaxby ay1nnxxh21()()()()()/2!nnnnny xy xhy xhy xh yx( )21(1)12!( )(, )(1)!ppnnnnnppnyyyyhyhhpyE x hhp公式中的導數(shù)用下面公式計算：公式中的導數(shù)用下面公式計算：例子：例子：2(,)(,)(,)(,)2(,)(,)(,)nnnnxnnynnnnxxnnxynnnyynnnynnnyf xyyfxyfxyyyfxyfxyyfxyyfxyy ( )cos04(0)1y xxxycos,sin,cosnnnnnny

22、xyxyx 231cossin/2cos/6nnnnnyyhxhxhx步長步長0.1 步長步長0.01龍格庫塔法：龍格庫塔法：對于常微分方程的數(shù)值解法，一個關(guān)鍵在于選擇精度高的算法對于常微分方程的數(shù)值解法，一個關(guān)鍵在于選擇精度高的算法計算下面公式中的積分。計算下面公式中的積分。要高精度的計算積分，常用的方法是適當添加計算節(jié)點，無妨要高精度的計算積分，常用的方法是適當添加計算節(jié)點，無妨設用設用 m 個節(jié)點計算上面積分，節(jié)點為個節(jié)點計算上面積分，節(jié)點為那么積分為那么積分為11()()( , ( )nnxnnxy xy xf x y x dx,1,2,inixxhim11( , ( )(, ()nn

23、mxininixif x y x dxhf xh y xh將積分改寫為：將積分改寫為：那么得公式：那么得公式：這樣的公式稱為顯式龍格庫塔公式。這樣的公式稱為顯式龍格庫塔公式。11( , ( )nnmxiixif x y x dxK11mnniiiyyK111(,)(,),2,3,nniininijjjKhf xyKhf xh yKim確定二階確定二階 RK 法法:系數(shù)滿足：系數(shù)滿足：此為四個未知數(shù)的三個方程，恣意取此為四個未知數(shù)的三個方程，恣意取 ，得，得11122122211(,)(,)nnnnnnyyKKKhf xyKhf xh yK12222121110,0,022 21222111,2

24、 特別，取特別，取 ，得到，得到通常也稱為變形歐拉法，也常寫為通常也稱為變形歐拉法，也常寫為它具有二階精度，也稱為二階二級它具有二階精度，也稱為二階二級 RK 方法。方法。1/212121(,)(/2,/2)nnnnnnyyKKhf xyKhf xhyK1(/2,(,)/2)nnnnnnyyhf xhyhf xy例子例子( )cos04(0)1y xxxy三階三級庫塔法三階三級庫塔法部分截斷誤差為部分截斷誤差為4階，整體截斷誤差為階，整體截斷誤差為3階階11231123121(4)6(,)(,)22(,2)nnnnnnnnyyKKKKhf xyKhKhf xyKhf xh yKK最常用的四級四

25、階公式，稱為龍格庫塔公式：最常用的四級四階公式，稱為龍格庫塔公式：部分截斷誤差為部分截斷誤差為5階，整體截斷誤差為階，整體截斷誤差為4階。階。1123411223431(22)6(,)(,)22(,)22(,)nnnnnnnnnnyyKKKKKhf xyKhKhf xyKhKhf xyKhf xh yK隱式龍格庫塔法：隱式龍格庫塔法：111221122(,)1,2,nnmmininiiimmyyKKKKhf xh yKKKim常用的有二階常用的有二階RK法：法：111111(,)22nnnnyyKKhf xh yK例子：例子：用隱式二階用隱式二階RK法：法：(01)(0)1yyxy 11111

26、0( , ),(/2)1/21/21nnnnnnhf x yyKh yKKyhhyyKyyhy 210(1/2)1nnyyhhy用顯式二階用顯式二階RK法：法： DOUBLE PRECISION h,y(0:100),yy(0:100)double precision ak OPEN(20,FILE=OUTPUT5.DAT,STATUS=UNKNOWN)h=1.0/100y(0)=1.0yy(0)=1.0do 10 i=1,100 y(i)=y(i-1)-(h/(1-h/2)*y(i-1) yy(i)=yy(i-1)*(1-h+h*h/2)10 continue do 20 i=0,100 w

27、rite(20,100) i*h,y(i),yy(i),exp(-i*h)100 format (1x,f4.2,f8.5,f8.5,f8.5)20 continueENDh=0.1h=0.01半隱式龍格庫塔法：半隱式龍格庫塔法：111221122,1111,11(,)(,)1,2,nnmmininiii iiiininii iiyyKKKKhf xh yKKKg K hfxh yKKim1122,11111,11(,)1(,)ininiii iiininii iiKhf xh yKKKhg fxh yKK最常用的二級三階半隱式龍格庫塔公式：最常用的二級三階半隱式龍格庫塔公式：11211121

28、10.413154321.4131543261 (1)(,)(,)661 (1)(,)6(,)0.17378667nnynnnnynnnnyyKKKhhfxyf xyKhhfxh yKf xh yK微分方程組微分方程組一階方程組一階方程組11122212121010202000( )( ,( ),( ),( )( )( ,( ),( ),( )( )( ,( ),( ),( )(),(),()mmmmmmmy xf x y xyxyxyxfx y xyxyxyxfx y xyxyxy xyyxyyxy00( )( , ( )()y xf x y xy xy 1212010200(,) ,(,)

29、 ,(,)mmmyy yyffffyyyy龍格庫塔公式龍格庫塔公式1123411223431(22)6(,)(,)22(,)22(,)nnnnnnnnnnyyKKKKKhf xyKhKhf xyKhKhf xyKhf xh yK寫成分量方式：寫成分量方式：,1,12341121112121221222312411322331(22)6(,)(,)2222(,)2222(,)1,2,i ni niiiiiinnnmnmiinnnmnmiinnnmniinnnmnmyyKKKKKhf xyyyKKKhKhf xyyyKKKhKhf xyyyKhf xh yKyKyKim例子：例子：1122212(

30、01)(0)0,(0)1yyyyyxyy 1,11,1112131411121121121212221312141132231(262)() ()22 ()22 ()nnnnnnnnnnyyKKKKKhyyKKKhyyKKKhyyKhyKyK2,12,212223242122122222232242231(262)()2()2()nnnnnnyyKKKKKhyKKh yKKh yKh yK DOUBLE PRECISION h,y1(0:100),y2(0:100)double precision ak1,ak2,ak3,ak4,bk1,bk2,bk3,bk4 OPEN(20,FILE=OUT

31、PUT6.DAT,STATUS=UNKNOWN)h=1.0/10y1(0)=0.0y2(0)=1.0do 10 i=1,10 ak1=h*(-y1(i-1)+y2(i-1) bk1=-h*y2(i-1) ak2=h*(-(y1(i-1)+ak1/2.0)+y2(i-1)+bk1/2.0) bk2=-h*(y2(i-1)+bk1/2.0) ak3=h*(-(y1(i-1)+ak2/2.0)+y2(i-1)+bk2/2.0) bk3=-h*(y2(i-1)+bk2/2.0) ak4=h*(-(y1(i-1)+ak3)+y2(i-1)+bk3) bk4=-h*(y2(i-1)+bk3) y1(i)=y1(i-1)+(ak1+2*ak2+2*ak3+ak4)/6.0 y2(i)=y2(i-1)+(bk1+2*bk2+2*bk3+bk4)/6.010 continue do 20 i=0,10 write(20,100) i*h,y1(i),i*h*exp(-i*h),y2(i),exp(-i*h)100 format (1x,f