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

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

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

3、 這個(gè)一階微分方程就不能用初等函數(shù)及其積分這個(gè)一階微分方程就不能用初等函數(shù)及其積分來(lái)表達(dá)它的解。來(lái)表達(dá)它的解。 引言引言3再如，方程再如，方程 1)0(yyy的解的解 , ,雖然有表可查雖然有表可查,但對(duì)于表上沒(méi)有給但對(duì)于表上沒(méi)有給出出 的值的值,仍需插值方法來(lái)計(jì)算仍需插值方法來(lái)計(jì)算xey xe引言引言4從實(shí)際問(wèn)題當(dāng)中歸納出來(lái)的微分方程，通常主要依從實(shí)際問(wèn)題當(dāng)中歸納出來(lái)的微分方程，通常主要依靠數(shù)值解法來(lái)解決。本章主要討論一階常微分方程靠數(shù)值解法來(lái)解決。本章主要討論一階常微分方程初值問(wèn)題初值問(wèn)題 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無(wú)關(guān)無(wú)關(guān)) )使使 引言引言5引言引言6對(duì)對(duì)R內(nèi)任意兩個(gè)內(nèi)任意兩個(gè) 都成立都成立, ,則方程則方程( 1 )的解的解 在在 a, b 上存在且唯一。上存在且唯一。 21, yy)(xyy 2121),(),(yyLyxfyxf9. 2 簡(jiǎn)單的數(shù)值方法簡(jiǎn)單的數(shù)值方法 對(duì)常微分方程初值問(wèn)題對(duì)常微分方程初值問(wèn)題( (1)式的數(shù)值解法，就是式的數(shù)值

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

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

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

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

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

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

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

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

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

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

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

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

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

18、xi-1)在在xi處展開(kāi)成處展開(kāi)成Taylor級(jí)數(shù)，即級(jí)數(shù)，即 )()(! 2)()()()(32111hOxxxyxxxyxyxyiiiiiiii 歐拉法的局部截?cái)嗾`差歐拉法的局部截?cái)嗾`差24)()(! 2)()(32hOxyhxyhxyiii ),(2)(11iiiiyxhfxyy由由 ),(2)()(! 2)()(32iiiiiyxhfhOxyhxyhxy )()(! 2)()(32hOxyhxyhxyiii 再將再將y(xi+1)在在xi+1處進(jìn)行泰勒展開(kāi)處進(jìn)行泰勒展開(kāi))()(! 2)()()(321hOxyhxyhxyxyiiii (8)(9)歐拉法的局部截?cái)嗾`差歐拉法的局部截?cái)嗾`差

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

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

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

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

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