第七章常微分方程數(shù)值解

1、計算方法河北大學(xué)電信學(xué)院第七章第七章 常微分方程的數(shù)值解法常微分方程的數(shù)值解法 7.1 7.1 引言引言 包含自變量、未知函數(shù)及未知函數(shù)的導(dǎo)數(shù)或微包含自變量、未知函數(shù)及未知函數(shù)的導(dǎo)數(shù)或微分的方程稱為微分方程。在微分方程中分的方程稱為微分方程。在微分方程中, , 自變量的自變量的個數(shù)只有一個個數(shù)只有一個, , 稱為常微分方程稱為常微分方程. .。自變量的個數(shù)。自變量的個數(shù)為兩個或兩個以上的微分方程叫偏微分方程。微分為兩個或兩個以上的微分方程叫偏微分方程。微分方程中出現(xiàn)的未知函數(shù)最高階導(dǎo)數(shù)的階數(shù)稱為微分方程中出現(xiàn)的未知函數(shù)最高階導(dǎo)數(shù)的階數(shù)稱為微分方程的階數(shù)。如果未知函數(shù)方程的階數(shù)。如果未知函數(shù)y

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

3、這個一階微分方程就不能用初等函數(shù)及其積這個一階微分方程就不能用初等函數(shù)及其積分來表達(dá)它的解。分來表達(dá)它的解。 再如，方程再如，方程 1)0(yyy的解的解 , ,雖然有表可查雖然有表可查, ,但對于表但對于表上沒有給出上沒有給出 的值的值, ,仍需插值方法來仍需插值方法來計算計算xey xe從實際問題當(dāng)中歸納出來的微分方程，通常主要依從實際問題當(dāng)中歸納出來的微分方程，通常主要依靠數(shù)值解法來解決。本章主要討論一階常微分方程靠數(shù)值解法來解決。本章主要討論一階常微分方程初值問題初值問題 00)(),(yxyyxfy ( 7.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滿足李普希茲滿足李普希茲( (lipschitzlipschitz) )條件，即存在常數(shù)條件，即存在常數(shù)l(l(它與它與x,yx,y無關(guān)無關(guān)) )使使 2121),(),(yylyxfyxf對對r內(nèi)任意兩個內(nèi)任意兩個 都成立都成立, ,則方程則方程( 7.1 )的解的解 在在 a, b 上存在且唯一。上存在且唯一。 21, yy)(xyy 7.2 數(shù)值方法的基本思想數(shù)值方法的基本思想 對常微分方程初值問題對常微分方程初值問題( (7.1)式的數(shù)值解法，就式的數(shù)值解法，就是要算出

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

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

7、行不同的離散化處理。 021,yyyyiii1iyy00)(),(yxyyxfy對于初值問題對于初值問題的數(shù)值解法，首先要解決的問題就是如何對微分方的數(shù)值解法，首先要解決的問題就是如何對微分方程進(jìn)行離散化，建立求數(shù)值解的遞推公式。遞推公程進(jìn)行離散化，建立求數(shù)值解的遞推公式。遞推公式通常有兩類，一類是計算式通常有兩類，一類是計算yi+1時只用到時只用到xi+1, xi 和和yi，即前一步的值，因此有了初值以后就可以逐步往下即前一步的值，因此有了初值以后就可以逐步往下計算，此類方法稱為單步法；其代表是計算，此類方法稱為單步法；其代表是龍格龍格庫塔庫塔法。另一類是計算法。另一類是計算y yi i+1

8、+1時，除用到時，除用到x xi i+1+1,x,xi i和和y yi i以外，以外，還要用到還要用到，即前面，即前面k步的值，此類方法稱為步的值，此類方法稱為多步法多步法；其代表；其代表是亞當(dāng)斯法。是亞當(dāng)斯法。 00)(),(yxyyxfy), 2 , 1( ,kpyxpipi7.3 歐拉（歐拉（euler）法）法7.3.1 euler公式公式 歐拉（歐拉（euler）方法是解初值問題的最）方法是解初值問題的最簡單的數(shù)值方法。初值問題簡單的數(shù)值方法。初值問題的解的解y=y(x)y=y(x)代表通過點代表通過點 的一條稱之為的一條稱之為微分方程的積分曲線。積分曲線上每一點微分方程的積分曲線。積

9、分曲線上每一點 的切線的斜率的切線的斜率 等于函數(shù)等于函數(shù) 在在這點的值。這點的值。 00)(),(yxyyxfy),(00yx),(yx)(xy),(yxf pi+1 pn y=y(x) p1 pi pn pi+1 p0 x0 x1 xi xi+1 xn pi p1 euler法的求解過程是法的求解過程是: :從初從初始點始點p0(即點即點( (x x0 0,y,y0 0)出發(fā)出發(fā), ,作積分曲線作積分曲線y=y(x)y=y(x)在在p0點上點上切線切線 ( (其斜率為其斜率為 ),),與與x=xx=x1 1直線直線10pp),()(000yxfxy相交于相交于p1點點( (即點即點( (x

10、 x1 1,y,y1 1),),得到得到y(tǒng) y1 1作為作為y(xy(x1 1) )的近似值的近似值, ,如上圖所示。過點如上圖所示。過點( (x x0 0,y,y0 0),),以以f(xf(x0 0,y,y0 0) )為為斜率的切線斜率的切線方程為方程為 當(dāng)當(dāng) 時時, ,得得 )(,(0000 xxyxfyy1xx )(,(010001xxyxfyy這樣就獲得了這樣就獲得了p p1 1點的坐標(biāo)。點的坐標(biāo)。 pi+1 pn y=y(x) p1 pi pn pi+1 p0 x0 x1 xi xi+1 xn pi p1 同樣同樣, 過過點點p1(x x1 1,y,y1 1),),作積分曲線作積分曲

11、線y=y(x)y=y(x)的切線的切線交直線交直線x=xx=x2 2于于p2點點, ,切線切線 的斜率的斜率 = =直線方程為直線方程為21pp)(1xy),(11yxf)(,(1111xxyxfyy)(,(121112xxyxfyy當(dāng)當(dāng) 時時, ,得得 2xx 當(dāng)當(dāng) 時時, ,得得 pi+1 pn y= y(x) p1 pi pn pi+1 p0 x0 x1 xi xi+1 xn pi p1 由此獲得了由此獲得了p p2 2的坐標(biāo)。重復(fù)以上過程的坐標(biāo)。重復(fù)以上過程, ,就可獲得一系就可獲得一系列的點列的點: :p p1 1, ,p p1 1, , ,p pn n。對已求得點對已求得點以以 =

12、 = 為斜率作直線為斜率作直線 ),(nnnyxp)(nxy),(nnyxf)(,(nnnnxxyxfyy1nxx)(,(11nxnnnnxxyxfyynnyxy)(取取 從圖形上看從圖形上看, ,就獲得了一條近似于曲線就獲得了一條近似于曲線y=y(x)y=y(x)的折線的折線 。 pi+ 1 pn y= y(x) p1 pi pn pi+ 1 p0 x0 x1 xi xi+ 1 xn pi p1 這樣這樣, ,從從x x0 0逐個算出逐個算出對應(yīng)的數(shù)值解對應(yīng)的數(shù)值解 nxxx,21nyyy,21npppp321通常取通常取 ( (常數(shù)常數(shù)),),則則euler法的計算格式法的計算格式 hhx

13、xiii1)(),(001xyyyxhfyyiiii i=0,1,n ( 7.2 ) 還可用數(shù)值微分、數(shù)值積分法和泰勒展開法推導(dǎo)還可用數(shù)值微分、數(shù)值積分法和泰勒展開法推導(dǎo)eulereuler格式。以數(shù)值積分為例進(jìn)行推導(dǎo)。格式。以數(shù)值積分為例進(jìn)行推導(dǎo)。將方程將方程 的兩端在區(qū)間的兩端在區(qū)間 上積分得，上積分得， ),(yxfy 1,iixx11),(iiiixxxxdxyxfdxy11)(,)(),()()(1iiiixxixxiidxxyxfxydxyxfxyxy 選擇不同的計算方法計算上式的積分項選擇不同的計算方法計算上式的積分項 , ,就會得到不同的計算公式。就會得到不同的計算公式。 1)

14、(,iixxdxxyxf(7.3) 用左矩形方法計算積分項用左矩形方法計算積分項 )(,)()(,11iiiixxxyxfxxdxxyxfii代入代入(7.3)(7.3)式式, ,并用并用y yi i近似代替式中近似代替式中y(xy(xi i) )即可得到即可得到向前歐拉（向前歐拉（eulereuler）公式）公式 ),(1iiiiyxhfyy 由于數(shù)值積分的矩形方法精度很低，所以由于數(shù)值積分的矩形方法精度很低，所以歐拉（歐拉（eulereuler）公式當(dāng)然很粗糙。）公式當(dāng)然很粗糙。 例例7.1 用歐拉法解初值問題用歐拉法解初值問題 1)0()6 . 00(2yxxyyy取步長取步長h=0.2

15、 ,h=0.2 ,計算過程保留計算過程保留4 4位小數(shù)位小數(shù) 解解: h=0.2, : h=0.2, 歐拉迭代格式歐拉迭代格式 2),(xyyyxf21),(iiiiiiiiyhxhyyyxhfyy)2 , 1 , 0()4(2 . 0iyxyiii當(dāng)當(dāng) k=0, x1=0.2時，已知時，已知x0=0，y0=1，有，有 y(0.2) y1=0.21(401)0.8當(dāng)當(dāng) k=1, x2=0.4時，已知時，已知x1 =0.2， y1 =0.8，有，有 y(0.4) y2 =0.20.8(40.20.8)0.6144當(dāng)當(dāng) k=2, x3 =0.6時，已知時，已知x2 =0.4， y2 =0.6144

16、，有，有 y(0.6) y3=0.20.6144(4-0.40.6144)=0.4613 7.3.2 梯形公式梯形公式為了提高精度為了提高精度, ,對方程對方程 的兩端在區(qū)間上的兩端在區(qū)間上 積分得，積分得，改用梯形方法計算其積分項，即改用梯形方法計算其積分項，即 ),(yxfy 1,iixx1)(,)()(1iixxiidxxyxfxyxy)(,()(,(2)(,1111iiiiiixxxyxfxyxfxxdxxyxfii( 7.4 ) 代入代入(7.4)(7.4)式式, ,并用近似代替式中即可得到梯形公式并用近似代替式中即可得到梯形公式 ),(),(2111iiiiiiyxfyxfhyy(

17、 7.5 ) 由于數(shù)值積分的梯形公式比矩形公式的精度高，由于數(shù)值積分的梯形公式比矩形公式的精度高，因此梯形公式（因此梯形公式（7.57.5）比歐拉公式）比歐拉公式( 7.2 )( 7.2 )的精度高的精度高一個數(shù)值方法。一個數(shù)值方法。 ),(),(2111iiiiiiyxfyxfhyy( 7.5 ) ( (7.5)式的右端含有未知的式的右端含有未知的y yi i+1+1, ,它是一個關(guān)于它是一個關(guān)于y yi i+1+1的函數(shù)方程的函數(shù)方程, ,這類數(shù)值方法稱為這類數(shù)值方法稱為隱式方法隱式方法。相。相反地反地, ,歐拉法是關(guān)于歐拉法是關(guān)于y yi i+1+1的一個直接的計算公式，的一個直接的計算

18、公式， 這類數(shù)值方法稱為這類數(shù)值方法稱為顯式方法。顯式方法。 7.3.3 兩步歐拉公式兩步歐拉公式 對方程對方程 的兩端在區(qū)間上的兩端在區(qū)間上 積分得積分得 ),(yxfy 11,iixx11)(,)()(11iixxiidxxyxfxyxy ( 7.6 ) 改用中矩形公式計算其積分項，即改用中矩形公式計算其積分項，即 )(,)(,1111iiiixxxyxfxxdxxyxfii代入上式代入上式, ,并用并用y yi i近似代替式中近似代替式中y(xy(xi i) )即可得到兩即可得到兩步歐拉公式步歐拉公式 ),(211iiiiyxhfyy ( 7.7 ) 前面介紹過的數(shù)值方法前面介紹過的數(shù)值

19、方法, ,無論是歐拉無論是歐拉方法方法, ,還是梯形方法，它們都是單步法還是梯形方法，它們都是單步法, ,其其特點是在計算特點是在計算y yi i+1+1時只用到前一步的信息時只用到前一步的信息y yi i; ;可是公式可是公式( (7.7)中除了中除了y yi i外外, ,還用到更前一步還用到更前一步的信息的信息y yi i-1-1, ,即調(diào)用了前兩步的信息即調(diào)用了前兩步的信息, ,故稱其故稱其為兩步歐拉公式為兩步歐拉公式 7.3.4 歐拉法的局部截斷誤差歐拉法的局部截斷誤差 衡量求解公式好壞的一個主要標(biāo)準(zhǔn)是求解公式的衡量求解公式好壞的一個主要標(biāo)準(zhǔn)是求解公式的精度精度, 因此引入局部截斷誤差

20、和階數(shù)的概念。因此引入局部截斷誤差和階數(shù)的概念。 定義定義7.1 在在yi準(zhǔn)確的前提下準(zhǔn)確的前提下, 即即 時時, 用數(shù)值用數(shù)值方法計算方法計算yi+1的誤差的誤差 , 稱為該數(shù)值方法稱為該數(shù)值方法計算時計算時yi+1的局部截斷誤差。的局部截斷誤差。 對于歐拉公式，假定對于歐拉公式，假定 ，則有，則有)(iixyy 11)(iiiyxyr)(iixyy )()()(,()(1iiiiiixyhxyxyxfhxyy而將真解而將真解y(x)在在xi處按二階泰勒展開處按二階泰勒展開 ),()(! 2)()()(121 iiiiixxyhxyhxyxy)(!2)(211yhyxyii 因此有因此有 定

21、義定義7.2 數(shù)值方法的局部截斷誤差為數(shù)值方法的局部截斷誤差為 ,則稱這則稱這種數(shù)值方法的階數(shù)是種數(shù)值方法的階數(shù)是p。步長。步長(h n 結(jié)束。結(jié)束。 10,xx)(21),(),(11cpipiiciiipyyyyxhfyyyxhfyy11, yx0101,yyxx（2）改進(jìn)歐拉法的流程圖）改進(jìn)歐拉法的流程圖 開 始 輸 入x0, y0,h , n 1 n x0 + h x1 y0+ h f( x0,y0 ) yp y0+ h f( x1,yp) yc ( yp+ yc) /2 y1 輸 出x1, y1 n + 1 n n = n ? x1 x0 y1 y0 結(jié) 束 n y ( (3)3)

22、程序?qū)崿F(xiàn)程序?qū)崿F(xiàn)( (見附錄見附錄a a-15 改進(jìn)歐拉法計算常微改進(jìn)歐拉法計算常微 分方程初值問題分方程初值問題 ) 例例7.2 7.2 用改進(jìn)歐拉法解初值問題用改進(jìn)歐拉法解初值問題 1)0(2yyxyy區(qū)間為區(qū)間為 0,10,1 , ,取步長取步長h=0.1h=0.1 解解: : 改進(jìn)歐拉法的具體形式改進(jìn)歐拉法的具體形式 )(21)2(1.0)2(1.011cpipipiciiiipyyyyxyyyyxyyy本題的精確解為本題的精確解為 , ,計算見計算見p p158158列表所示列表所示 xxy21)(例例7.3 對初值問題對初值問題 1)0(0yyy證明用梯形公式求得的近似解為證明用梯

23、形公式求得的近似解為 nnhhynhx 并證明當(dāng)步長并證明當(dāng)步長h h0 0時時, ,y yn n收斂于精確解收斂于精確解證明證明: : 解初值問題的梯形公式為解初值問題的梯形公式為 xe),(),(nnnnnnyxfyxfhyyyyxf),( 211nnnnyyhyy 整理成顯式整理成顯式 nnyhhy反復(fù)迭代反復(fù)迭代, ,得到得到 yhhyhhyhhyhhynnnnn.10ynnhhy22 由于由于 ，有，有 nhx xxxxhxhhhxhnhhhhhyeee2121lim22limlim222222000nnhhy22xnhy elim0 證畢證畢 7.4 7.4 龍格龍格- -庫塔（庫

24、塔（runge-kuttarunge-kutta）法）法7.4.1 7.4.1 龍格龍格- -庫塔庫塔( (runge-kuttarunge-kutta) )法的基本思想法的基本思想 eulereuler公式可改寫成公式可改寫成 ),(111iiiiyxfkhkyy則則yi i+1+1的表達(dá)式的表達(dá)式y(tǒng)( (xi i+1+1) )與的與的taylortaylor展開式的前兩項展開式的前兩項完全相同完全相同, ,即局部截斷誤差為即局部截斷誤差為 。 改進(jìn)的改進(jìn)的eulereuler公式又可改寫成公式又可改寫成 )(2ho),(),()(21121211hkyxfkyxfkkkhyyiiiiii

25、上述兩組公式在形式上有一個共同點上述兩組公式在形式上有一個共同點: :都是用都是用f(x,y)f(x,y)在某些點上值的線性組合得出在某些點上值的線性組合得出y(xy(xi i+1+1) )的近似的近似值值y yi i+1+1, ,而且增加計算的次數(shù)而且增加計算的次數(shù)f f( (x x, ,y y) )的次數(shù)的次數(shù), ,可提高截可提高截斷誤差的階。如歐拉公式斷誤差的階。如歐拉公式: :每步計算一次每步計算一次f f( (x x, ,y y) )的的值值, ,為一階方法。改進(jìn)歐拉公式需計算兩次為一階方法。改進(jìn)歐拉公式需計算兩次f f( (x x, ,y y) )的值，它是二階方法。它的局部截斷誤

26、差為的值，它是二階方法。它的局部截斷誤差為 。)(3ho 于是可考慮用函數(shù)于是可考慮用函數(shù)f( (x, ,y) )在若干點上的函在若干點上的函數(shù)值的線性組合來構(gòu)造近似公式，構(gòu)造時要求數(shù)值的線性組合來構(gòu)造近似公式，構(gòu)造時要求近似公式在近似公式在( (xi i, ,yi i) )處的處的taylor展開式與解展開式與解y( (x) )在在xi i處的處的taylor展開式的前面幾項重合，從而展開式的前面幾項重合，從而使近似公式達(dá)到所需要的階數(shù)。既避免求偏導(dǎo)使近似公式達(dá)到所需要的階數(shù)。既避免求偏導(dǎo), ,又提高了計算方法精度的階數(shù)?；蛘哒f又提高了計算方法精度的階數(shù)?；蛘哒f, ,在在 這一步內(nèi)多預(yù)報幾個

27、點的斜率值，然后這一步內(nèi)多預(yù)報幾個點的斜率值，然后將其加權(quán)平均作為平均斜率，則可構(gòu)造出更高將其加權(quán)平均作為平均斜率，則可構(gòu)造出更高精度的計算格式，這就是龍格精度的計算格式，這就是龍格庫塔（庫塔（runge-kutta）法的基本思想。）法的基本思想。 1,iixx7.4.2 二階龍格二階龍格庫塔法庫塔法 在在 上取兩點上取兩點xi i和和 , ,以該兩點處的以該兩點處的斜率值斜率值k1 1和和k2 2的加權(quán)平均的加權(quán)平均( (或稱為線性組合或稱為線性組合) )來求取平來求取平均斜率均斜率k* *的近似值的近似值k，即，即 1,iixxphxxipi2211kkk式中式中: :k1 1為為xi i

28、點處的切線斜率值，點處的切線斜率值， k2 2為為 點處的切線斜率值點處的切線斜率值, ,比照改進(jìn)的歐比照改進(jìn)的歐拉法拉法, ,將將 視為視為 ，即可得，即可得 )(),(1iiixyyxfkpixphxi1ix),(12phkyphxfkii對常微分方程初值問題對常微分方程初值問題(7.1)(7.1)式的解式的解 y= =y( (x),),根據(jù)微根據(jù)微分中值定理，存在點分中值定理，存在點 ，使得，使得 ),(1iixx)(,()(yfyk式中式中 k k可看作是可看作是y=y(x)y=y(x)在區(qū)間在區(qū)間 上的平均斜率。上的平均斜率。所以可得計算公式為：所以可得計算公式為： 1,iixxhk

29、xyxyii)()(1)()(2211kkhxyi（7.14） 將將y(xy(xi i) )在在x=xx=xi i處進(jìn)行二階處進(jìn)行二階taylortaylor展開：展開： )()(! 2)()()(321hoxyhxyhxyxyiiii （7.15） )()()(11iiiixxyxyxy也即也即 hkxyxyii)()(1（7.13）將將 在在x=xx=xi i處進(jìn)行一階處進(jìn)行一階taylortaylor展開：展開： ),()(12phkyphxfphxykiii)(),(),(),(),(22hoyxfyxfyxfphyxfkiiyiiiixii)()()(2hoxyphxyii 將以上結(jié)

30、果代入（將以上結(jié)果代入（7.147.14）得：）得： )()()(22111kkhxyxyii)()()()()(221hoxyphxyxyhxyiiii )()()()()(32221hoxyphxyhxyiii （7.16） 對式對式(7.15)(7.15)和和(7.16)(7.16)進(jìn)行比較系數(shù)后可知進(jìn)行比較系數(shù)后可知, ,只要只要 211221p（7.17） 成立成立, ,格式格式(7.14)(7.14)的局部截斷誤差就等于的局部截斷誤差就等于)(3ho有有2 2階階精度精度式式(7.17)(7.17)中具有三個未知量中具有三個未知量, ,但只有兩個方程但只有兩個方程, ,因而因而有無

31、窮多解。若取有無窮多解。若取 , ,則則p p=1=1，這是無窮多解，這是無窮多解中的一個解，將以上所解的值代入式中的一個解，將以上所解的值代入式(7.14)(7.14)并改寫并改寫可得可得 2121),(),()(21121211hkyxfkyxfkkkhyyiiiiii 不難發(fā)現(xiàn)，上面的格式就是改進(jìn)的歐拉格式。不難發(fā)現(xiàn)，上面的格式就是改進(jìn)的歐拉格式。凡滿足條件式（凡滿足條件式（7.177.17）有一簇形如上式的計算格式，）有一簇形如上式的計算格式，這些格式統(tǒng)稱為二階龍格這些格式統(tǒng)稱為二階龍格庫塔格式。因此改進(jìn)的庫塔格式。因此改進(jìn)的歐拉格式是眾多的二階龍格歐拉格式是眾多的二階龍格庫塔法中的一

32、種特殊庫塔法中的一種特殊格式。格式。 若取若取 , ,則則 ，此時二階龍格，此時二階龍格- -庫塔庫塔法的計算公式為法的計算公式為 0121, 12p)2,(),(1212121khyxfkyxfkhkyyiiiiii1,2 , 1 , 0ni 此計算公式稱為變形的二階龍格此計算公式稱為變形的二階龍格庫塔法。式中庫塔法。式中 為區(qū)間為區(qū)間 的中點。的中點。 21ix1,iixx7.4.3 三階龍格三階龍格- -庫塔法庫塔法 為了進(jìn)一步提高精度，設(shè)除為了進(jìn)一步提高精度，設(shè)除 外再增加一點外再增加一點 pix) 1(qpqhxxiqi并用三個點并用三個點 , , , , 的斜率的斜率k1 1, ,

33、k2 2, ,k3 3加權(quán)平均加權(quán)平均得出平均斜率得出平均斜率k* *的近似值，這時計算格式具有形式的近似值，這時計算格式具有形式: : ixpixqix),(),()1 (1213211phkyphxfkyxfkkkkhyyiiiiii（7.18） 為了預(yù)報點為了預(yù)報點 的斜率值的斜率值k3 3, ,在區(qū)間在區(qū)間 內(nèi)有兩內(nèi)有兩個斜率值個斜率值k1 1和和k2 2可以用可以用, ,可將可將k1 1, ,k2 2加權(quán)平均得出加權(quán)平均得出 上的平均斜率上的平均斜率, ,從而得到從而得到 的預(yù)報值的預(yù)報值 qixqiixx,qiixx,)(qixyqiy21)1 (kkqhyyiqi于是可得于是可得

34、 ),(3qiqiyxfk運用運用taylor展開方法選擇參數(shù)展開方法選擇參數(shù) , ,可以使可以使格式格式( (7.18)的局部截斷誤差為的局部截斷誤差為 , ,即具有三階精度，即具有三階精度，這類格式統(tǒng)稱為這類格式統(tǒng)稱為三階龍格三階龍格庫塔方法庫塔方法。下列是其中。下列是其中的一種，稱為的一種，稱為庫塔（庫塔（kutta）公式）公式。 ,qp)(4ho)4(6)2(,()2,(),(3211211312121kkkhyykkhyxfkkhyxfkyxfkiiiiiiii（7.19） 7.4.4 四階龍格四階龍格庫塔法庫塔法 如果需要再提高精度，用類似上述的處理方法，如果需要再提高精度，用類似

35、上述的處理方法，只需在區(qū)間只需在區(qū)間 上用四個點處的斜率加權(quán)平均作為上用四個點處的斜率加權(quán)平均作為平均斜率平均斜率k*的近似值，構(gòu)成一系列四階龍格的近似值，構(gòu)成一系列四階龍格庫塔公庫塔公式。具有四階精度，即局部截斷誤差是式。具有四階精度，即局部截斷誤差是 。 由于推導(dǎo)復(fù)雜，這里從略，只介紹最常用的一種由于推導(dǎo)復(fù)雜，這里從略，只介紹最常用的一種四階經(jīng)典龍格四階經(jīng)典龍格庫塔公式庫塔公式。 qiixx,)(5ho)22(6),()2,()2,(),(43211314221312121kkkkhyyhkyxfkkhyxfkkhyxfkyxfkiiiiiiiiii（7.20） 7.4.5 7.4.5 四

36、階龍格四階龍格庫塔法算法實現(xiàn)庫塔法算法實現(xiàn)(1)(1) 計算步驟計算步驟 輸入輸入 , ,h,nh,n 使用龍格使用龍格庫塔公式（庫塔公式（7.207.20）計算出）計算出y y1 1 輸出輸出 ，并使，并使 轉(zhuǎn)到轉(zhuǎn)到 直至直至n n n n 結(jié)束。結(jié)束。 10,xx11, yx0101,yyxx(2) (2) 四階龍格四階龍格庫塔算法流程圖庫塔算法流程圖 開 始 輸 入 x0, y0,h , n 1 n x0 + h x1 f(x0,y0 ) k1, f(x0+ h /2 ,y0 + h k1/2 ) k2 f(x0+ h /2 ,y0 + h k2/2 ) k3, f(x1,y0+ h k

37、3) k4 y0+ h (k1+ 2 k2+ + 2 k3+ k4)/6 y1 輸 出 x1, y1 n + 1 n n = n ? x1 x0 y1 y0 結(jié) 束 n y (3)(3) 程序?qū)崿F(xiàn)程序?qū)崿F(xiàn)( (見附錄見附錄a a-16 四階龍格四階龍格- -庫塔法計庫塔法計 算常微分方程初值問題算常微分方程初值問題) ) 例例7.4 7.4 取步長取步長h=0.2h=0.2，用經(jīng)典格式求解初值問題，用經(jīng)典格式求解初值問題 1)0(2yxyy10 x解解: : 由四階龍格由四階龍格- -庫塔公式可得庫塔公式可得 2 . 0, 1, 0,2),(00hyxxyyxf0),(001yxfk2 . 0

38、) 1 , 1 . 0()2,(10202fkhyxfkh204. 0)02. 1 , 1 . 0()2,(20203fkhyxfkh41632. 0)0408. 1 , 2 . 0(),(3004fhkyxfkh040811. 1)22(62 . 0432101kkkkyy可同樣進(jìn)行其余可同樣進(jìn)行其余y yi i的計算。本例方程的解為的計算。本例方程的解為，數(shù)值解，數(shù)值解y yi i與準(zhǔn)確解與準(zhǔn)確解y(xy(xi i) )的對照表見教材的對照表見教材p p163163所示，所示，從表中看到所求的數(shù)值解具有從表中看到所求的數(shù)值解具有4位有效數(shù)字。位有效數(shù)字。 2xey 龍格龍格庫塔方法的推導(dǎo)基

39、于庫塔方法的推導(dǎo)基于taylor展開方法，展開方法，因而它要求所求的解具有較好的光滑性。如果解的因而它要求所求的解具有較好的光滑性。如果解的光滑性差，那么，使用四階龍格光滑性差，那么，使用四階龍格庫塔方法求得的庫塔方法求得的數(shù)值解，其精度可能反而不如改進(jìn)的歐拉方法。在數(shù)值解，其精度可能反而不如改進(jìn)的歐拉方法。在實際計算時，應(yīng)當(dāng)針對問題的具體特點選擇合適的實際計算時，應(yīng)當(dāng)針對問題的具體特點選擇合適的算法。算法。 7.4.6 變步長的龍格變步長的龍格-庫塔法庫塔法 在微分方程的數(shù)值解中，選擇適當(dāng)?shù)牟介L是非常在微分方程的數(shù)值解中，選擇適當(dāng)?shù)牟介L是非常重要的。單從每一步看，步長越小，截斷誤差就越重要的

40、。單從每一步看，步長越小，截斷誤差就越??；但隨著步長的縮小，在一定的求解區(qū)間內(nèi)所要?。坏S著步長的縮小，在一定的求解區(qū)間內(nèi)所要完成的步數(shù)就增加了。這樣會引起計算量的增大，完成的步數(shù)就增加了。這樣會引起計算量的增大，并且會引起舍入誤差的大量積累與傳播。因此微分并且會引起舍入誤差的大量積累與傳播。因此微分方程數(shù)值解法也有選擇步長的問題。方程數(shù)值解法也有選擇步長的問題。 以經(jīng)典的四階龍格以經(jīng)典的四階龍格- -庫塔法庫塔法( (7.20)為例。從節(jié)點為例。從節(jié)點x xi i出發(fā)，先以出發(fā)，先以h為步長求出一個近似值，記為為步長求出一個近似值，記為 ，由于局部截斷誤差為由于局部截斷誤差為 ，故有，故有

41、)(1hiy)(5ho5)(11)(chyxyhii當(dāng)h h值不大時，式中的系數(shù)值不大時，式中的系數(shù)c c可近似地看作為常數(shù)?？山频乜醋鳛槌?shù)。然后將步長折半然后將步長折半, ,即以為即以為 步長步長, ,從節(jié)點從節(jié)點x xi i出發(fā)出發(fā), ,跨跨兩步到節(jié)點兩步到節(jié)點x xi i+1+1, ,再求得一個近似值再求得一個近似值 , ,每跨一步的每跨一步的截斷誤差是截斷誤差是 , ,因此有因此有2h)2(1hiy52 hc5)2(1122)()(hcxyxyhii這樣這樣 161)()()(11)2(11hiihiiyxyyxy)(151)()(1)2(1)2(11hihihiiyyyxy由此可

42、得由此可得 這表明以這表明以 作為作為 的近似值，其誤差可用步的近似值，其誤差可用步長折半前后兩次計算結(jié)果的偏差長折半前后兩次計算結(jié)果的偏差 )2(1hiy)(1ixy)(1)2(1hihiyy來判斷所選步長是否適當(dāng)來判斷所選步長是否適當(dāng)當(dāng)要求的數(shù)值精度為當(dāng)要求的數(shù)值精度為時：時： （1 1）如果）如果，反復(fù)將步長折半進(jìn)行計算，直，反復(fù)將步長折半進(jìn)行計算，直至至為止為止, ,并取其最后一次步長的計算結(jié)果作為并取其最后一次步長的計算結(jié)果作為 （2 2）如果）如果為止，并以上一次步長的計算結(jié)果作為為止，并以上一次步長的計算結(jié)果作為 。 這種通過步長加倍或折半來處理步長的方法稱為這種通過步長加倍或折

43、半來處理步長的方法稱為變步長法。表面上看，為了選擇步長，每一步都要變步長法。表面上看，為了選擇步長，每一步都要反復(fù)判斷反復(fù)判斷，增加了計算工作量，但在方程的解，增加了計算工作量，但在方程的解y(x)y(x)變化劇烈的情況下，總的計算工作量得到減少，結(jié)變化劇烈的情況下，總的計算工作量得到減少，結(jié)果還是合算的。果還是合算的。1iy1iy7.5 亞當(dāng)姆斯方法亞當(dāng)姆斯方法7.5.1 亞當(dāng)姆斯格式亞當(dāng)姆斯格式 龍格龍格-庫塔方法是一類重要算法，但這類庫塔方法是一類重要算法，但這類算法在每一步都需要先預(yù)報幾個點上的斜率算法在每一步都需要先預(yù)報幾個點上的斜率值，計算量比較大。考慮到計算值，計算量比較大?？紤]

44、到計算yi+1之前已得之前已得出一系列節(jié)點上出一系列節(jié)點上 的斜率值，的斜率值，能否利用這些已知值來減少計算量呢？能否利用這些已知值來減少計算量呢？ 這就是亞當(dāng)姆斯（這就是亞當(dāng)姆斯（adams）方法的設(shè)計）方法的設(shè)計思想。思想。 11,xxxii 設(shè)用設(shè)用x xi i,x,xi i+1+1兩點的斜率值加權(quán)平均作為區(qū)間兩點的斜率值加權(quán)平均作為區(qū)間 上的平均斜率，有計算格式上的平均斜率，有計算格式 1,iixx),(),()1 (11111iiiiiiiiiiyxfyyxfyyyhyy（7.21） 選取參數(shù)選取參數(shù)，使格式（，使格式（7.217.21）具有二階精度。）具有二階精度。 將將 在在x

45、xi i處處taylor展開展開 1iy)()(! 21)(321hohyhyyyiii 代入計算格式代入計算格式(7.21)(7.21)化簡化簡, ,并假設(shè)并假設(shè), , 因此有因此有 )(),(11iiiixyyxyy)()()()(321hoxyhxyhxyyiiii 與與y(xi+1)在在xi處的處的taylor展開式展開式 )()(21)()()(321hoxyhxyhxyxyiiii 相比較相比較, 需取需取 21才使格式才使格式(7.21)具有二階精度。這樣導(dǎo)出的計算格式具有二階精度。這樣導(dǎo)出的計算格式 )3(211iiiiyyhyy稱之為二階亞當(dāng)姆斯格式。類似地可以導(dǎo)出三階亞稱之

46、為二階亞當(dāng)姆斯格式。類似地可以導(dǎo)出三階亞當(dāng)姆斯格式當(dāng)姆斯格式。 )51623(12211iiiiiyyyhyy和和四階亞當(dāng)姆斯格式。四階亞當(dāng)姆斯格式。 )9375955(243211iiiiiiyyyyhyy( (7.22) 這里和下面均記這里和下面均記 ),(kikikiyxfy 上述幾種亞當(dāng)姆斯格式都是顯式的，算法比較上述幾種亞當(dāng)姆斯格式都是顯式的，算法比較簡單，但用節(jié)點簡單，但用節(jié)點 的斜率值來預(yù)報的斜率值來預(yù)報區(qū)間區(qū)間 上的平均斜率是個外推過程，效果不夠上的平均斜率是個外推過程，效果不夠理想。為了進(jìn)一步改善精度，變外推為內(nèi)插，即理想。為了進(jìn)一步改善精度，變外推為內(nèi)插，即增加節(jié)點增加節(jié)點

47、xi+1的斜率值來得出的斜率值來得出 上的平均斜率。上的平均斜率。譬如考察形如譬如考察形如 11,xxxii1,iixx1,iixx),(),()1 (11111iiiiiiiiiiyxfyyxfyyyhyy（7.23） 的隱式格式的隱式格式, ,設(shè)設(shè)( (7.23)右端的右端的taylor展開有展開有 )(),(11iiiixyyxyy)()()1 ()()(321hoxyhxyhxyyiiii 可見要使格式可見要使格式(7.23)(7.23)具有二階精度具有二階精度, ,需令需令 , ,這樣就可構(gòu)造二階隱式亞當(dāng)姆斯格式這樣就可構(gòu)造二階隱式亞當(dāng)姆斯格式 21)(211iiiiyyhyy其實是

48、梯形格式。類似可導(dǎo)出三階隱式亞當(dāng)姆斯格其實是梯形格式。類似可導(dǎo)出三階隱式亞當(dāng)姆斯格式式 )85(12111iiiiiyyyhyy和四階隱式亞當(dāng)姆斯格式和四階隱式亞當(dāng)姆斯格式 )5199(242111iiiiiiyyyyhyy（7.247.24） 7.5.2 亞當(dāng)姆斯預(yù)報亞當(dāng)姆斯預(yù)報-校正格式校正格式 參照改進(jìn)的歐拉格式的構(gòu)造方法，以四階亞當(dāng)參照改進(jìn)的歐拉格式的構(gòu)造方法，以四階亞當(dāng)姆斯為例，將顯式（姆斯為例，將顯式（7.22）和隱式（）和隱式（7.24）相結(jié)合，）相結(jié)合，用顯式公式做預(yù)報，再用隱式公式做校正，可構(gòu)成用顯式公式做預(yù)報，再用隱式公式做校正，可構(gòu)成亞當(dāng)姆斯預(yù)報亞當(dāng)姆斯預(yù)報-校正格式校正

49、格式 )9375955(243211iiiiiiyyyyhyy),(111iiiyxfy)5199 (242111iiiiiiyyyyhyy),(111iiiyxfy（7.257.25） 預(yù)報：預(yù)報： 校正：校正： 這種預(yù)報這種預(yù)報- -校正格式是四步法，它在計校正格式是四步法，它在計算算y yi i+1+1時不但用到前一步的信息時不但用到前一步的信息 ，而且，而且要用到再前面三步的信息要用到再前面三步的信息 ，因，因此它不能自行啟動。在實際計算時，可借助此它不能自行啟動。在實際計算時，可借助于某種單步法，譬如四階龍格于某種單步法，譬如四階龍格庫塔法提供庫塔法提供開始值開始值 。 iiyy,3

50、21,iiiyyy321,yyy例例7.5 取步長取步長h=0.1,h=0.1,用亞當(dāng)姆斯預(yù)報用亞當(dāng)姆斯預(yù)報- -校正公式求解校正公式求解 初值問題初值問題 1)0(yyxy10 x的數(shù)值解。的數(shù)值解。 解解: : 用四階龍格用四階龍格- -庫塔公式求出發(fā)值庫塔公式求出發(fā)值 , ,計算得：計算得：1 . 0, 1, 0,),(00hyxyxyxf321,yyy399717. 1,242805. 1,110342. 1321yyy表中的表中的 , ,y yi i和和y(xy(xi i) )分別為預(yù)報值、校正值和準(zhǔn)確分別為預(yù)報值、校正值和準(zhǔn)確解解( ),( ),以比較計算結(jié)果的精度。以比較計算結(jié)果

51、的精度。 再使用亞當(dāng)姆斯預(yù)報再使用亞當(dāng)姆斯預(yù)報- -校正公式校正公式(7.25),(7.25),見教材見教材p p166166列表算得其余的計算結(jié)果列表算得其余的計算結(jié)果iy12xeyx 7.6 算法的穩(wěn)定性及收斂性算法的穩(wěn)定性及收斂性 7.6.1 穩(wěn)定性穩(wěn)定性 穩(wěn)定性在微分方程的數(shù)值解法中是一個非常穩(wěn)定性在微分方程的數(shù)值解法中是一個非常重要的問題。因為微分方程初值問題的數(shù)值方法是重要的問題。因為微分方程初值問題的數(shù)值方法是用差分格式進(jìn)行計算的，而在差分方程的求解過程用差分格式進(jìn)行計算的，而在差分方程的求解過程中，存在著各種計算誤差，這些計算誤差如舍入誤中，存在著各種計算誤差，這些計算誤差如舍

52、入誤差等引起的擾動，在傳播過程中，可能會大量積累，差等引起的擾動，在傳播過程中，可能會大量積累，對計算結(jié)果的準(zhǔn)確性將產(chǎn)生影響。這就涉及到算法對計算結(jié)果的準(zhǔn)確性將產(chǎn)生影響。這就涉及到算法穩(wěn)定性問題。穩(wěn)定性問題。 當(dāng)在某節(jié)點上當(dāng)在某節(jié)點上x xi i的的y yi i值有大小為值有大小為的擾動時，的擾動時，如果在其后的各節(jié)點如果在其后的各節(jié)點 上的值上的值y yi i產(chǎn)生的偏差產(chǎn)生的偏差都不大于都不大于，則稱這種方法是穩(wěn)定的。，則稱這種方法是穩(wěn)定的。 穩(wěn)定性不僅與算法有關(guān)，而且與方程中函數(shù)穩(wěn)定性不僅與算法有關(guān)，而且與方程中函數(shù)f(x,y)f(x,y)也有關(guān)，討論起來比較復(fù)雜。也有關(guān)，討論起來比較復(fù)雜

53、。為簡單起見，為簡單起見，通常只針對模型方程通常只針對模型方程)(ijxj)0(yy來討論。一般方程若局部線性化，也可化為上述形來討論。一般方程若局部線性化，也可化為上述形式。模型方程相對比較簡單，若一個數(shù)值方法對模式。模型方程相對比較簡單，若一個數(shù)值方法對模型方程是穩(wěn)定的，并不能保證該方法對任何方程都型方程是穩(wěn)定的，并不能保證該方法對任何方程都穩(wěn)定，但若某方法對模型方程都不穩(wěn)定，也就很難穩(wěn)定，但若某方法對模型方程都不穩(wěn)定，也就很難用于其他方程的求解。用于其他方程的求解。先考察顯式先考察顯式eulereuler方法的穩(wěn)定性。模型方程方法的穩(wěn)定性。模型方程的的eulereuler公式為公式為 y

54、y)(),(1iiiiiiyhyyxhfyy), 2 , 1 , 0()1 (iyhi將上式反復(fù)遞推后，可得將上式反復(fù)遞推后，可得 011)1(yhyii), 2 , 1()1 (00iyyhyiii或或h1式中式中 要使要使y yi i有界，其充要條件為有界，其充要條件為 111h即即 由于由于 ，故有，故有 020 h（7.267.26） 可見，如欲保證算法的穩(wěn)定，顯式可見，如欲保證算法的穩(wěn)定，顯式eulereuler格式格式的步長的步長h h的選取要受到式（的選取要受到式（7.267.26）的限制。）的限制。 的絕的絕對值越大，則限制的對值越大，則限制的h h值就越小。值就越小。 用隱式

55、用隱式eulereuler格式，對模型方程格式，對模型方程 的計算公式為，可化為的計算公式為，可化為)(11iiiyhyyiiyhy111由于由于 , ,則恒有則恒有 , ,故恒有故恒有 0111hiiyy1 因此，隱式因此，隱式eulereuler格式是絕對穩(wěn)定的（無條件穩(wěn)格式是絕對穩(wěn)定的（無條件穩(wěn)定）的（對任何定）的（對任何h0h0）。）。 7.6.2 7.6.2 收斂性收斂性 常微分方程初值問題的求解，是將微分方程轉(zhuǎn)化常微分方程初值問題的求解，是將微分方程轉(zhuǎn)化為差分方程來求解，并用計算值為差分方程來求解，并用計算值y yi i來近似替代來近似替代y(xy(xi i) )，這種近似替代是否

56、合理，還須看分割區(qū)間這種近似替代是否合理，還須看分割區(qū)間 的的長度長度h h越來越小時，即越來越小時，即 時，時， 是是否成立。若成立，則稱該方法是收斂的，否則稱為否成立。若成立，則稱該方法是收斂的，否則稱為不收斂。不收斂。 這里仍以這里仍以eulereuler方法為例，來分析其收斂性。方法為例，來分析其收斂性。eulereuler格式格式iixx,101iixxh)(iixyy),(1iiiiyxhfyy設(shè)設(shè) 表示取表示取 時時, 按按euler公式的計算結(jié)果公式的計算結(jié)果, 即即1iy)(iixyy )(,()(1iiiixyxhfxyyeuler方法局部截斷誤差為：方法局部截斷誤差為：)

57、()(2)(1211 iiiixxyhyxy設(shè)有常數(shù)設(shè)有常數(shù) ，則，則 )(max21xycbxa 211)(chyxyii（7.27） 總體截斷誤差總體截斷誤差 1111111)()(iiiiiiiyyyxyyxy)(,()(),(11iiiiiiiixyxhfxyyxhfyyy),()(,()(iiiiiiyxfxyxfhyxy（7.28） 又又 由于由于f(x,y)f(x,y)關(guān)于關(guān)于y y滿足李普希茨條件滿足李普希茨條件, ,即即 iiiiiiyxylyxfxyxf)(),()(代入代入(7.28)上式，有上式，有 iiiiyxyhlyy)()1(11ihl)1( 再利用式（再利用式（

58、7.27），式（），式（7.28） 1111111)()(iiiiiiiyyyxyyxy2)1 (chhli21)1 (chhlii即即 上式反復(fù)遞推后，可得上式反復(fù)遞推后，可得 1020)1 ()1 (ikkiihlchhl1)1 ()1 (0iihllchhl（7.29） 設(shè)設(shè) （t t為常數(shù)）為常數(shù)） tihxxi0hlehl 1tlihlieehl)1 (因為因為 所以所以 把上式代入式（把上式代入式（7.297.29），得），得 ) 1(0tltlielche若不計初值誤差，即若不計初值誤差，即 ，則有，則有 00helctli)1(（7.30） 式式(7.30)(7.30)說明說明

59、, ,當(dāng)時當(dāng)時 , , ,即即 , ,所以所以eulereuler方法是收斂的，且其收斂速度為方法是收斂的，且其收斂速度為 ，即具，即具有一階收斂速度。同時還說明有一階收斂速度。同時還說明eulereuler方法的整體截斷方法的整體截斷誤差為誤差為 ，因此算法的精度為一階。，因此算法的精度為一階。 0h0i)(iixyy )(ho)(ho7.7 一階常微分方程組與高階方程一階常微分方程組與高階方程 我們已介紹了一階常微分方程初值問題的各種我們已介紹了一階常微分方程初值問題的各種數(shù)值解法，對于一階常微分方程組，可類似得到各數(shù)值解法，對于一階常微分方程組，可類似得到各種解法，而高階常微分方程可轉(zhuǎn)化

60、為一階常微分方種解法，而高階常微分方程可轉(zhuǎn)化為一階常微分方程組來求解。程組來求解。 7.7.1 一階常微分方程組一階常微分方程組對于一階常微分方程組的初值問題對于一階常微分方程組的初值問題 0000)(),()(),(zxzzyxgzyxyzyxfy（7.31） 可以把單個方程可以把單個方程 中的中的f 和和y看作向量看作向量來處理，這樣就可把前面介紹的各種差分算法推廣來處理，這樣就可把前面介紹的各種差分算法推廣到求一階方程組初值問題中來。到求一階方程組初值問題中來。 ),(yxfy 設(shè)設(shè) 為節(jié)點上的近似解，為節(jié)點上的近似解，則有改進(jìn)的則有改進(jìn)的eulereuler格式為格式為 iiizyii