常微分方程初值問題數(shù)值解法

1、第第9章章 常微分方程初值問題數(shù)值解法常微分方程初值問題數(shù)值解法9.1 引言引言9.2 簡單的數(shù)值方法簡單的數(shù)值方法9.1引言引言 包含自變量、未知函數(shù)及未知函數(shù)的導數(shù)或微包含自變量、未知函數(shù)及未知函數(shù)的導數(shù)或微分的方程稱為微分方程。在微分方程中分的方程稱為微分方程。在微分方程中, 自變量的自變量的個數(shù)只有一個個數(shù)只有一個, 稱為稱為常微分方程常微分方程。自變量的個數(shù)為。自變量的個數(shù)為兩個或兩個以上的微分方程叫兩個或兩個以上的微分方程叫偏微分方程偏微分方程。微分方。微分方程中出現(xiàn)的未知函數(shù)最高階導數(shù)的階數(shù)稱為微分方程中出現(xiàn)的未知函數(shù)最高階導數(shù)的階數(shù)稱為微分方程的階數(shù)。如果未知函數(shù)程的階數(shù)。如果

2、未知函數(shù)y及其各階導數(shù)及其各階導數(shù)都是一次的都是一次的,則稱它是線性的則稱它是線性的,否則稱為非線性的。否則稱為非線性的。 )(,nyyy 2 在高等數(shù)學中，對于常微分方程的求解，給出在高等數(shù)學中，對于常微分方程的求解，給出了一些典型方程求解析解的基本方法，如了一些典型方程求解析解的基本方法，如可分離變可分離變量法量法、常系數(shù)齊次線性方程常系數(shù)齊次線性方程的解法、常系數(shù)非齊次的解法、常系數(shù)非齊次線性方程的解法等。但能求解的常微分方程仍然是線性方程的解法等。但能求解的常微分方程仍然是有限的，大多數(shù)的常微分方程是不可能給出解析解。有限的，大多數(shù)的常微分方程是不可能給出解析解。 譬如譬如 22yxy

3、 這個一階微分方程就不能用初等函數(shù)及其積分這個一階微分方程就不能用初等函數(shù)及其積分來表達它的解。來表達它的解。 引言引言3再如，方程再如，方程 1)0(yyy的解的解 , ,雖然有表可查雖然有表可查,但對于表上沒有給但對于表上沒有給出出 的值的值,仍需插值方法來計算仍需插值方法來計算xey xe引言引言4從實際問題當中歸納出來的微分方程，通常主要依從實際問題當中歸納出來的微分方程，通常主要依靠數(shù)值解法來解決。本章主要討論一階常微分方程靠數(shù)值解法來解決。本章主要討論一階常微分方程初值問題初值問題 00)(),(yxyyxfy ( 1 ) 在區(qū)間在區(qū)間a x ba x b上的數(shù)值解法上的數(shù)值解法。

4、 可以證明可以證明, ,如果函數(shù)在帶形區(qū)域如果函數(shù)在帶形區(qū)域 R=axb,R=axb,-y y 內(nèi)連續(xù)，且關(guān)于內(nèi)連續(xù)，且關(guān)于y y滿足李普希茲滿足李普希茲(Lipschitz)(Lipschitz)條件條件，即存在常數(shù)，即存在常數(shù)L(L(它與它與x,yx,y無關(guān)無關(guān)) )使使 引言引言5引言引言6對對R內(nèi)任意兩個內(nèi)任意兩個 都成立都成立, ,則方程則方程( 1 )的解的解 在在 a, b 上存在且唯一。上存在且唯一。 21, yy)(xyy 2121),(),(yyLyxfyxf9. 2 簡單的數(shù)值方法簡單的數(shù)值方法 對常微分方程初值問題對常微分方程初值問題( (1)式的數(shù)值解法，就是式的數(shù)值

5、解法，就是要算出精確解要算出精確解y(x)y(x)在區(qū)間在區(qū)間 a,b 上的一系列離散節(jié)點上的一系列離散節(jié)點 處的函數(shù)值處的函數(shù)值 的近似值的近似值 。相鄰兩個節(jié)點的間距相鄰兩個節(jié)點的間距 稱為步長，步長可稱為步長，步長可以相等，也可以不等。本章總是假定以相等，也可以不等。本章總是假定h為定數(shù)，稱為定數(shù)，稱為為定步長定步長，這時節(jié)點可表示為，這時節(jié)點可表示為bxxxxann110)(,),(),(10nxyxyxynyyy,10iixxh1niihxxi, 2 , 1 ,07數(shù)值方法的基本思想數(shù)值方法的基本思想8數(shù)值解法需要把連續(xù)性的問題加以離散化，從數(shù)值解法需要把連續(xù)性的問題加以離散化，從而

6、求出離散節(jié)點的數(shù)值解。而求出離散節(jié)點的數(shù)值解。 對常微分方程數(shù)值解法的基本出發(fā)點就是離散對常微分方程數(shù)值解法的基本出發(fā)點就是離散化。其數(shù)值解法有兩個基本特點，它們都采用化。其數(shù)值解法有兩個基本特點，它們都采用“步步進式進式”，即求解過程順著節(jié)點排列的次序一步一步，即求解過程順著節(jié)點排列的次序一步一步地向前推進，描述這類算法，要求給出用已知信息地向前推進，描述這類算法，要求給出用已知信息 計算計算 的遞推公式。建立這的遞推公式。建立這類遞推公式的基本方法是在這些節(jié)點上用數(shù)值積分、類遞推公式的基本方法是在這些節(jié)點上用數(shù)值積分、021,yyyyiii1iy中的導數(shù)中的導數(shù) 進行不同的離散化處理進行不

7、同的離散化處理。對于初值問題對于初值問題y00)(),(yxyyxfy數(shù)值方法的基本思想數(shù)值方法的基本思想900)(),(yxyyxfy數(shù)值微分、泰勒展開等離散化方法，對初值問題數(shù)值微分、泰勒展開等離散化方法，對初值問題的數(shù)值解法，首先要解決的問題就是的數(shù)值解法，首先要解決的問題就是如何對微分方如何對微分方程進行離散化程進行離散化，建立求數(shù)值解的，建立求數(shù)值解的遞推公式遞推公式。遞推公。遞推公式通常有兩類，一類是計算式通常有兩類，一類是計算yi+1時只用到時只用到xi+1, xi 和和yi，即前一步的值，因此有了初值以后就可以逐步往下即前一步的值，因此有了初值以后就可以逐步往下計算，此類方法稱

8、為計算，此類方法稱為單步法單步法；其代表是；其代表是龍格龍格庫塔庫塔法。另一類是計算法。另一類是計算y yi+1i+1時，除用到時，除用到x xi+1i+1,x,xi i和和y yi i以外，以外，還要用到還要用到 ，即前面，即前面k步的值，此類步的值，此類方法稱為方法稱為多步法多步法；其代表是；其代表是亞當斯法亞當斯法。 ), 2 , 1( ,kpyxpipi數(shù)值方法的基本思想數(shù)值方法的基本思想101 Euler公式公式 歐拉（歐拉（Euler）方法是解初值問題的最簡單的）方法是解初值問題的最簡單的數(shù)值方法。初值問題數(shù)值方法。初值問題的解的解y=y(x)代表通過點代表通過點 的一條稱之為的一

9、條稱之為微分方微分方程的積分曲線程的積分曲線。積分曲線上每一點。積分曲線上每一點 的切線的切線的斜率的斜率 等于函數(shù)等于函數(shù) 在這點的值。在這點的值。 00)(),(yxyyxfy),(00yx),(yx)(xy),(yxf11Euler公式公式12通常取通常取 ( (常數(shù)常數(shù)),),hhxxiii1將方程將方程 的兩端在區(qū)間的兩端在區(qū)間 上積分得，上積分得， ),(yxfy 1,iixx11),(iiiixxxxdxyxfdxy11)(,)(),()()(1iiiixxixxiidxxyxfxydxyxfxyxy(1) 選擇不同的計算方法計算上式的積分項選擇不同的計算方法計算上式的積分項 ,

10、就會得到不同的計算公式。就會得到不同的計算公式。 1)(,iixxdxxyxfEuler公式公式13用用左矩形方法左矩形方法計算積分項計算積分項 )(,)()(,11iiiixxxyxfxxdxxyxfii 代入代入(1)式式,并用并用yi近似代替式中近似代替式中y(xi)即可得到向前歐拉即可得到向前歐拉（Euler）公式）公式 ),(1iiiiyxhfyy 由于數(shù)值積分的矩形方法精度很低，所以歐拉由于數(shù)值積分的矩形方法精度很低，所以歐拉（Euler）公式當然很粗糙。）公式當然很粗糙。 Euler公式公式14Euler公式公式15Euler法的計算格式法的計算格式 )(),(001xyyyxh

11、fyyiiii i=0,1,n ( 2 ) 例例1 用歐拉法解初值問題用歐拉法解初值問題 1)0()6 . 00(2yxxyyy取步長取步長h=0.2 ,計算過程保留計算過程保留4位小數(shù)位小數(shù) 解解: h=0.2, 歐拉迭代格式歐拉迭代格式 2),(xyyyxf21),(iiiiiiiiyhxhyyyxhfyy)2 , 1 , 0()4(2 . 0iyxyiii例題例題16例題例題17當當 i=0, x1=0.2時，已知時，已知x0=0，y0=1，有，有 y(0.2) y1=0.21(401)0.8當當 i=1, x2=0.4時，已知時，已知x1 =0.2， y1 =0.8，有，有 y(0.4

12、) y2 =0.20.8(40.20.8)0.6144當當 i=2, x3 =0.6時，已知時，已知x2 =0.4， y2 =0.6144，有，有 y(0.6) y3=0.20.6144(4-0.40.6144)=0.4613 2 梯形公式梯形公式 為了提高精度為了提高精度, ,對方程對方程 的兩端在區(qū)間上的兩端在區(qū)間上 積分得，積分得， 改用梯形方法計算其積分項，即改用梯形方法計算其積分項，即 ),(yxfy 1,iixx1)(,)()(1iixxiidxxyxfxyxy)(,()(,(2)(,1111iiiiiixxxyxfxyxfxxdxxyxfii( 4 ) 代入代入(4)(4)式式,

13、 ,并用近似代替式中即可得到梯形公式并用近似代替式中即可得到梯形公式 梯形公式梯形公式18 (5)式的右端含有未知的式的右端含有未知的yi+1,它是一個關(guān)于它是一個關(guān)于yi+1的函數(shù)方程的函數(shù)方程,這類數(shù)值方法稱為這類數(shù)值方法稱為隱式方法隱式方法。相反地。相反地,歐拉法是關(guān)于歐拉法是關(guān)于yi+1的一個直接的計算公式，的一個直接的計算公式， 這類數(shù)這類數(shù)值方法稱為顯式方法。值方法稱為顯式方法。 梯形公式梯形公式19 由于數(shù)值積分的梯形公式比矩形公式的精度高，由于數(shù)值積分的梯形公式比矩形公式的精度高，因此梯形公式（因此梯形公式（5）比歐拉公式）比歐拉公式(2)的精度高一個數(shù)值的精度高一個數(shù)值方法。

14、方法。 ),(),(2111iiiiiiyxfyxfhyy(5 ) 3 兩步歐拉公式兩步歐拉公式 對方程對方程 的兩端在區(qū)間上的兩端在區(qū)間上 積分得積分得 ),(yxfy 11,iixx11)(,)()(11iixxiidxxyxfxyxy ( 6 ) 改用中矩形公式計算其積分項，即改用中矩形公式計算其積分項，即 )(,)(,1111iiiixxxyxfxxdxxyxfii代入上式代入上式,并用并用yi近似代替式中近似代替式中y(xi)即可得到兩步歐即可得到兩步歐拉公式拉公式 兩步歐拉公式兩步歐拉公式20 前面介紹過的數(shù)值方法前面介紹過的數(shù)值方法,無論是歐拉方法無論是歐拉方法,還還是梯形方法，

15、它們都是單步法是梯形方法，它們都是單步法,其特點是在計算其特點是在計算yi+1時只用到前一步的信息時只用到前一步的信息yi;可是公式可是公式(7)中除了中除了yi外外,還用到更前一步的信息還用到更前一步的信息yi-1,即調(diào)用了前兩步的信息即調(diào)用了前兩步的信息,故稱其為兩步歐拉公式故稱其為兩步歐拉公式 兩步歐拉公式兩步歐拉公式21),(211iiiiyxhfyy (7 ) 4 歐拉法的局部截斷誤差歐拉法的局部截斷誤差 衡量求解公式好壞的一個主要標準是求解公式的衡量求解公式好壞的一個主要標準是求解公式的精度精度, 因此引入局部截斷誤差和階數(shù)的概念。因此引入局部截斷誤差和階數(shù)的概念。 定義定義1 在

16、在yi準確的前提下準確的前提下, 即即 時時, 用數(shù)值用數(shù)值方法計算方法計算yi+1的誤差的誤差 , 稱為該數(shù)值方稱為該數(shù)值方法計算時法計算時yi+1的局部截斷誤差。的局部截斷誤差。 對于歐拉公式，假定對于歐拉公式，假定 ，則有，則有)(iixyy 11)(iiiyxyR)(iixyy )()()(,()(1iiiiiixyhxyxyxfhxyy歐拉法的局部截斷誤差歐拉法的局部截斷誤差22歐拉法的局部截斷誤差歐拉法的局部截斷誤差23而將真解而將真解y(x)在在xi處按二階泰勒展開處按二階泰勒展開 ),()(! 2)()()(121 iiiiixxyhxyhxyxy因此有因此有 )(!2)(21

17、1yhyxyii 定義定義2 數(shù)值方法的局部截斷誤差為數(shù)值方法的局部截斷誤差為 ,則稱這種則稱這種數(shù)值方法的階數(shù)是數(shù)值方法的階數(shù)是P。步長。步長(h1) 越小越小,P越高越高, 則局則局部截斷誤差越小部截斷誤差越小,計算精度越高。歐拉公式的局部截計算精度越高。歐拉公式的局部截斷誤差為斷誤差為 , 歐拉方法僅為一階方法。歐拉方法僅為一階方法。)(1phO)(2hO兩步歐拉公式比歐拉公式精度也是高一個數(shù)值方法兩步歐拉公式比歐拉公式精度也是高一個數(shù)值方法,設(shè)設(shè) , 前兩步準確前兩步準確,則兩步歐拉則兩步歐拉公式公式 )(iixyy )(11iixyy)(,(2)(11iiiixyxhfxyy把把y(

18、xi-1)在在xi處展開成處展開成Taylor級數(shù)，即級數(shù)，即 )()(! 2)()()()(32111hOxxxyxxxyxyxyiiiiiiii 歐拉法的局部截斷誤差歐拉法的局部截斷誤差24)()(! 2)()(32hOxyhxyhxyiii ),(2)(11iiiiyxhfxyy由由 ),(2)()(! 2)()(32iiiiiyxhfhOxyhxyhxy )()(! 2)()(32hOxyhxyhxyiii 再將再將y(xi+1)在在xi+1處進行泰勒展開處進行泰勒展開)()(! 2)()()(321hOxyhxyhxyxyiiii (8)(9)歐拉法的局部截斷誤差歐拉法的局部截斷誤差

19、25歐拉法的局部截斷誤差歐拉法的局部截斷誤差26所以所以, 由由(8)和和(9)可得兩步歐拉公式的局部截斷誤差可得兩步歐拉公式的局部截斷誤差為為, 即即 )()(311hOyxyii5 改進的歐拉公式改進的歐拉公式 顯式歐拉公式計算工作量小顯式歐拉公式計算工作量小,但精度低。梯形公但精度低。梯形公式雖提高了精度式雖提高了精度,但為隱式公式但為隱式公式,需用迭代法求解需用迭代法求解,計算計算工作量大。綜合歐拉公式和梯形公式便可得到改進工作量大。綜合歐拉公式和梯形公式便可得到改進的歐拉公式。的歐拉公式。 先用歐拉公式先用歐拉公式(2)求出一個初步的近似值求出一個初步的近似值 ,稱稱為預(yù)測值為預(yù)測值

20、, 它的精度不高它的精度不高, 再用梯形公式再用梯形公式(5)對它對它1iy改進的歐拉公式改進的歐拉公式27改進的歐拉公式改進的歐拉公式28),(),(2),(1111iiiiiiiiiiyxfyxfhyyyxhfyy預(yù)測預(yù)測 校正校正 （10） 1,2 , 1 , 0ni 可以證明可以證明,公式公式(10)的精度為二階。這是一種一的精度為二階。這是一種一步顯式格式步顯式格式,它可以表示為嵌套形式它可以表示為嵌套形式。校正一次校正一次,即迭代一次即迭代一次,求得求得yi+1,稱為校正值稱為校正值, 這種預(yù)這種預(yù)測測-校正方法稱為改進的歐拉公式：校正方法稱為改進的歐拉公式：),(,(),(211

21、iiiiiiiiyxhfyxfyxfhyy（11）或者表示成下列平均化形式或者表示成下列平均化形式 )(21),(),(11cpipiiciiipyyyyxhfyyyxhfyy（12） 1,2 , 1 , 0ni改進的歐拉公式改進的歐拉公式29例例2 用改進歐拉法解初值問題用改進歐拉法解初值問題 1)0(2yyxyy區(qū)間為區(qū)間為 0,1 ,取步長取步長h=0.1 解解: 改進歐拉法的具體形式改進歐拉法的具體形式 )(21)2(1.0)2(1.011cpipipiciiiipyyyyxyyyyxyyy本題的精確解為本題的精確解為 ,計算結(jié)果略。計算結(jié)果略。 xxy21)(例題例題309. 3 龍

22、格龍格-庫塔方法庫塔方法1 龍格龍格-庫塔庫塔(Runge-Kutta)法法的基本思想的基本思想 Euler公式可改寫成公式可改寫成 ),(111iiiiyxfKhKyy則則yi+1的表達式的表達式y(tǒng)(xi+1)與的與的Taylor展開式的前兩項完展開式的前兩項完全相同全相同,即局部截斷誤差為即局部截斷誤差為 。 改進的改進的Euler公式又可改寫成公式又可改寫成 )(2hO31 上述兩組公式在形式上有一個共同點上述兩組公式在形式上有一個共同點:都是用都是用f(x,y)在某些點上值的線性組合得出在某些點上值的線性組合得出y(xi+1)的近似值的近似值yi+1,而且增加計算的次數(shù)而且增加計算的次

23、數(shù)f(x,y)的次數(shù)的次數(shù),可提高截斷誤可提高截斷誤差的階。如歐拉公式差的階。如歐拉公式:每步計算一次每步計算一次f(x,y)的值的值,為一為一階方法。改進歐拉公式需計算兩次階方法。改進歐拉公式需計算兩次f(x,y)的值，它是的值，它是二階方法。它的局部截斷誤差為二階方法。它的局部截斷誤差為 。)(3hO龍格龍格-庫塔法的基本思想庫塔法的基本思想32),(),()(21121211hKyxfKyxfKKKhyyiiiiii 于是可考慮用函數(shù)于是可考慮用函數(shù)f(x,y)在若干點上的函數(shù)值的在若干點上的函數(shù)值的線性組合來構(gòu)造線性組合來構(gòu)造近似公式近似公式，構(gòu)造時要求近似公式在，構(gòu)造時要求近似公式在(xi,yi)處的處的Taylor展開式與解展開式與解y(x)在在xi處的處的Taylor展開式展開式的前面幾項重合，從而使近似公式達到所需要的階數(shù)。的前面幾項重合，從而使近似公式達到所需要的階數(shù)。既避免求偏導既避免求偏導,又提高了計算方法精度的階數(shù)?；蛘哂痔岣吡擞嬎惴椒ň鹊碾A數(shù)?；蛘哒f說,在在 這一步內(nèi)多預(yù)報幾個點的斜率值，然后這一步內(nèi)多預(yù)報幾個點的斜率值，然后將其加權(quán)平均作為平均斜率，則可構(gòu)造出更高精