第9微分方程初值問題數值解法

1、第第9章章 常微分方程初值問題數值解法常微分方程初值問題數值解法9.1 9.1 引言引言 包含自變量、未知函數及未知函數的導數或微包含自變量、未知函數及未知函數的導數或微分的方程稱為微分方程。在微分方程中分的方程稱為微分方程。在微分方程中, , 自變量的自變量的個數只有一個個數只有一個, , 稱為常微分方程稱為常微分方程. .。自變量的個數。自變量的個數為兩個或兩個以上的微分方程叫偏微分方程。微分為兩個或兩個以上的微分方程叫偏微分方程。微分方程中出現的未知函數最高階導數的階數稱為微分方程中出現的未知函數最高階導數的階數稱為微分方程的階數。如果未知函數方程的階數。如果未知函數y y及其各階導數及

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

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

4、域如果函數在帶形區(qū)域 R=axb,R=axb,-y y 內連續(xù)，且關于內連續(xù)，且關于y y滿足李普希茲滿足李普希茲(Lipschitz)(Lipschitz)條件，即存在常數條件，即存在常數L(L(它與它與x,yx,y無關無關) )使使 2121),(),(yyLyxfyxf對對R內任意兩個內任意兩個 都成立都成立, ,則方程則方程( 9.1 )的解的解 在在 a, b 上存在且唯一。上存在且唯一。 21, yy)(xyy 數值方法的基本思想數值方法的基本思想 對常微分方程初值問題對常微分方程初值問題( (9.1)式的數值解法，就式的數值解法，就是要算出精確解是要算出精確解y(x)y(x)在區(qū)

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

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

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

8、,x,xi i和和y yi i以外，以外，還要用到還要用到，即前面，即前面k步的值，此類方法稱為步的值，此類方法稱為多步法多步法；其代表；其代表是亞當斯法。是亞當斯法。 00)(),(yxyyxfy), 2 , 1( ,kpyxpipi9.2 簡單的數值方法與基本概念簡單的數值方法與基本概念9.2.1 Euler公式公式 歐拉（歐拉（Euler）方法是解初值問題的最）方法是解初值問題的最簡單的數值方法。初值問題簡單的數值方法。初值問題的解的解y=y(x)y=y(x)代表通過點代表通過點 的一條稱之為的一條稱之為微分方程的積分曲線。積分曲線上每一點微分方程的積分曲線。積分曲線上每一點 的切線的斜

9、率的切線的斜率 等于函數等于函數 在在這點的值。這點的值。 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(x1 1,y,y1 1),),得

10、到得到y(tǒng) y1 1作為作為y(xy(x1 1) )的近似值的近似值, ,如上圖所示。過點如上圖所示。過點(x(x0 0,y,y0 0),),以以f(xf(x0 0,y,y0 0) )為為斜率的切線斜率的切線方程為方程為 當當 時時, ,得得 )(,(0000 xxyxfyy1xx )(,(010001xxyxfyy這樣就獲得了這樣就獲得了P P1 1點的坐標。點的坐標。 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),),作積分曲線作積分曲線y=y(x)y=y(x)的切線的切線

11、交直線交直線x=xx=x2 2于于P2點點, ,切線切線 的斜率的斜率 = =直線方程為直線方程為21PP)(1xy),(11yxf)(,(1111xxyxfyy)(,(121112xxyxfyy當當 時時, ,得得 2xx 當當 時時, ,得得 Pi+1 Pn y= y(x) P1 Pi Pn Pi+1 P0 x0 x1 xi xi+1 xn Pi P1 由此獲得了由此獲得了P P2 2的坐標。重復以上過程的坐標。重復以上過程, ,就可獲得一系就可獲得一系列的點列的點: :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逐個算出逐個算出對應的數值解對應的數值解 nxxx,21nyyy,21nPPPP321通常取通常取 ( (常數常數),),則則Euler法的計算格式法的計算格式 hhxxiii1)(),(001xyyyxh

13、fyyiiii i=0,1,n ( 9.2 ) 還可用還可用數值微分數值微分、數值積分法數值積分法和和泰勒展開法泰勒展開法推導推導EulerEuler格式。以數值積分為例進行推導。格式。以數值積分為例進行推導。將方程將方程 的兩端在區(qū)間的兩端在區(qū)間 上積分得，上積分得， ),(yxfy 1,iixx11),(iiiixxxxdxyxfdxy11)(,)(),()()(1iiiixxixxiidxxyxfxydxyxfxyxy 選擇不同的計算方法計算上式的積分項選擇不同的計算方法計算上式的積分項 , ,就會得到不同的計算公式。就會得到不同的計算公式。 1)(,iixxdxxyxf(9.3) 用左

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

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

16、6144(4-0.40.6144)=0.4613 0.1,2, 01) (0)1, hEulerxyyxyy 取利用公式求解（例例 . .2 . 01.1)2( 1計算及結果如下解：nnnnnnnnyxyyxyhyyclear; y=1, x=0, %初始化for n=1:10 y=1.1*y-0.2*x/y, x=x+0.1,endy = 1 x = 0 y = 1.1000 x = 0.1000y = 1.1918 x = 0.2000y = 1.2774 x = 0.3000y = 1.3582 x = 0.4000y = 1.4351 x = 0.5000y = 1.5090 x =

17、0.6000y = 1.5803 x = 0.7000y = 1.6498 x = 0.8000y = 1.7178 x = 0.9000y = 1.7848 x = 1.00009.2.2 梯形公式梯形公式為了提高精度為了提高精度, ,對方程對方程 的兩端在區(qū)間上的兩端在區(qū)間上 積分得，積分得，改用梯形方法計算其積分項，即改用梯形方法計算其積分項，即 ),(yxfy 1,iixx1)(,)()(1iixxiidxxyxfxyxy)(,()(,(2)(,1111iiiiiixxxyxfxyxfxxdxxyxfii(9.4 ) 代入代入(7.4)(7.4)式式, ,并用近似代替式中即可得到梯形公

18、式并用近似代替式中即可得到梯形公式 ),(),(2111iiiiiiyxfyxfhyy( 9.5 ) 由于數值積分的梯形公式比矩形公式的精度高，由于數值積分的梯形公式比矩形公式的精度高，因此梯形公式（因此梯形公式（9.59.5）比歐拉公式）比歐拉公式( 9.2 )( 9.2 )的精度高的精度高一個數值方法。一個數值方法。 ),(),(2111iiiiiiyxfyxfhyy( 9.5 ) ( (9.5)式的右端含有未知的式的右端含有未知的y yi+1i+1, ,它是一個關于它是一個關于y yi+1i+1的函數方程的函數方程, ,這類數值方法稱為這類數值方法稱為隱式方法隱式方法。相。相反地反地,

19、,歐拉法是關于歐拉法是關于y yi+1i+1的一個直接的計算公式，的一個直接的計算公式， 這類數值方法稱為這類數值方法稱為顯式方法。顯式方法。 9.2.3 兩步歐拉公式兩步歐拉公式 對方程對方程 的兩端在區(qū)間上的兩端在區(qū)間上 積分得積分得 ),(yxfy 11,iixx11)(,)()(11iixxiidxxyxfxyxy ( 9.6 ) 改用中矩形公式計算其積分項，即改用中矩形公式計算其積分項，即 )(,)(,1111iiiixxxyxfxxdxxyxfii代入上式代入上式, ,并用并用y yi i近似代替式中近似代替式中y(xy(xi i) )即可得到兩即可得到兩步歐拉公式步歐拉公式 ),

20、(211iiiiyxhfyy ( 9.7 ) 前面介紹過的數值方法前面介紹過的數值方法, ,無論是歐拉無論是歐拉方法方法, ,還是梯形方法，它們都是單步法還是梯形方法，它們都是單步法, ,其其特點是在計算特點是在計算y yi+1i+1時只用到前一步的信息時只用到前一步的信息y yi i; ;可是公式可是公式( (7.7)中除了中除了y yi i外外, ,還用到更前一步還用到更前一步的信息的信息y yi-1i-1, ,即調用了前兩步的信息即調用了前兩步的信息, ,故稱其故稱其為兩步歐拉公式為兩步歐拉公式 9.2.4 歐拉法的局部截斷誤差歐拉法的局部截斷誤差 衡量求解公式好壞的一個主要標準是求解公

21、式的衡量求解公式好壞的一個主要標準是求解公式的精度精度, 因此引入局部截斷誤差和階數的概念。因此引入局部截斷誤差和階數的概念。 定義定義9.1 在在yi準確準確的前提下的前提下, 即即 時時, 用數值用數值方法計算方法計算yi+1的誤差的誤差 , 稱為該數值方法稱為該數值方法計算時計算時yi+1的局部截斷誤差。的局部截斷誤差。 對于歐拉公式，假定對于歐拉公式，假定 ，則有，則有)(iixyy 11)(iiiyxyR)(iixyy )()()(,()(1iiiiiixyhxyxyxfhxyy而將真解而將真解y(x)在在xi處按二階泰勒展開處按二階泰勒展開 ),()(! 2)()()(121 ii

22、iiixxyhxyhxyxy)(!2)(211yhyxyii 因此有因此有 定義定義9.2 數值方法的局部截斷誤差為數值方法的局部截斷誤差為 ,則稱這則稱這種數值方法的階數是種數值方法的階數是P。步長步長(h N 結束。結束。 10,xx)(21),(),(11cpipiiciiipyyyyxhfyyyxhfyy11, yx0101,yyxx（2）改進歐拉法的流程圖）改進歐拉法的流程圖 開 始 輸 入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

23、n n = N ? x1 x0 y1 y0 結 束 n y ( (3)3) 程序實現程序實現( (改進歐拉法計算常微改進歐拉法計算常微 分方程初值問題分方程初值問題 ) 例例9.2 9.2 用改進歐拉法解初值問題用改進歐拉法解初值問題 1)0(2yyxyy區(qū)間為區(qū)間為 0,10,1 , ,取步長取步長h=0.1h=0.1 解解: : 改進歐拉法的具體形式改進歐拉法的具體形式 )(21)2(1.0)2(1.011cpipipiciiiipyyyyxyyyyxyyy本題的精確解為本題的精確解為 , ,xxy21)(clearx=0,yn=1 %初始化for n=1:10yp=yn+0.1*(yn-

24、2*x/yn); %預測x=x+0.1;yc=yn+0.1*(yp-2*x/yp) ;yn=(yp+yc)/2 %校正end例例9.3 對初值問題對初值問題 1)0(0yyy證明用梯形公式求得的近似解為證明用梯形公式求得的近似解為 nnhhynhx 并證明當步長并證明當步長h h0 0時時,y,yn n收斂于精確解收斂于精確解證明證明: : 解初值問題的梯形公式為解初值問題的梯形公式為 xe),(),(nnnnnnyxfyxfhyyyyxf),( 211nnnnyyhyy 整理成顯式整理成顯式 nnyhhy反復迭代反復迭代, ,得到得到 yhhyhhyhhyhhynnnnn.10ynnhhy2

25、2 由于由于 ，有，有 nhx xxxxhxhhhxhnhhhhhyeee2121lim22limlim222222000nnhhy22xnhy elim0 證畢證畢 9.3 9.3 龍格龍格- -庫塔（庫塔（Runge-KuttaRunge-Kutta）法）法9.3.1 9.3.1 龍格龍格- -庫塔庫塔(Runge-Kutta)(Runge-Kutta)法的基本思想法的基本思想 Euler Euler公式可改寫成公式可改寫成 ),(111iiiiyxfKhKyy則則yi+1i+1的表達式的表達式y(tǒng)( (xi+1i+1) )與的與的TaylorTaylor展開式的前兩項展開式的前兩項完全相同

26、完全相同, ,即局部截斷誤差為即局部截斷誤差為 。 改進的改進的EulerEuler公式又可改寫成公式又可改寫成 )(2hO),(),()(21121211hKyxfKyxfKKKhyyiiiiii 上述兩組公式在形式上有一個共同點上述兩組公式在形式上有一個共同點: :都是用都是用f(x,y)f(x,y)在某些點上值的線性組合得出在某些點上值的線性組合得出y(xy(xi+1i+1) )的近似的近似值值y yi+1i+1, ,而且增加計算的次數而且增加計算的次數f f( (x x, ,y y) )的次數的次數, ,可提高截可提高截斷誤差的階。如歐拉公式斷誤差的階。如歐拉公式: :每步計算一次每步

27、計算一次f f( (x x, ,y y) )的的值值, ,為一階方法。改進歐拉公式需計算兩次為一階方法。改進歐拉公式需計算兩次f f( (x x, ,y y) )的值，它是二階方法。它的局部截斷誤差為的值，它是二階方法。它的局部截斷誤差為 。)(3hO 于是可考慮用函數于是可考慮用函數f( (x, ,y) )在若干點上的函在若干點上的函數值的線性組合來構造近似公式，構造時要求數值的線性組合來構造近似公式，構造時要求近似公式在近似公式在( (xi i, ,yi i) )處的處的Taylor展開式與解展開式與解y( (x) )在在xi i處的處的Taylor展開式的前面幾項重合，從而展開式的前面幾

28、項重合，從而使近似公式達到所需要的階數。既避免求偏導使近似公式達到所需要的階數。既避免求偏導, ,又提高了計算方法精度的階數?；蛘哒f又提高了計算方法精度的階數?；蛘哒f, ,在在 這一步內多這一步內多預報幾個點的斜率值，然后預報幾個點的斜率值，然后將其加權平均作為平均斜率將其加權平均作為平均斜率，則可構造出更高，則可構造出更高精度的計算格式，這就是龍格精度的計算格式，這就是龍格庫塔（庫塔（Runge-Kutta）法的基本思想。）法的基本思想。 1,iixx9.3.2 二階龍格二階龍格庫塔法庫塔法 在在 上取兩點上取兩點xi i和和 , ,以該兩點處的以該兩點處的斜率值斜率值k1 1和和k2 2的

29、加權平均的加權平均( (或稱為線性組合或稱為線性組合) )來求取平來求取平均斜率均斜率k* *的近似值的近似值K，即，即 1,iixxphxxipi2211kkK式中式中: :k1 1為為xi i點處的切線斜率值，點處的切線斜率值， k2 2為為 點處的切線斜率值點處的切線斜率值, ,比照改進的歐比照改進的歐拉法拉法, ,將將 視為視為 ，即可得，即可得 )(),(1iiixyyxfkpixphxi1ix),(12phkyphxfkii對常微分方程初值問題對常微分方程初值問題(9.1)(9.1)式的解式的解 y= =y( (x),),根據微根據微分中值定理，存在點分中值定理，存在點 ，使得，使

30、得 ),(1iixx)(,()(yfyK式中式中 K K可看作是可看作是y=y(x)y=y(x)在區(qū)間在區(qū)間 上的平均斜率。上的平均斜率。所以可得計算公式為：所以可得計算公式為： 1,iixxhKxyxyii)()(1)()(2211kkhxyi（9.14） 將將y(xy(xi i) )在在x=xx=xi i處進行二階處進行二階TaylorTaylor展開：展開： )()(! 2)()()(321hOxyhxyhxyxyiiii （9 9.15） )()()(11iiiixxyxyxy也即也即 hKxyxyii)()(1（9.13）將將 在在x=xx=xi i處進行一階處進行一階TaylorT

31、aylor展開：展開： ),()(12phkyphxfphxykiii)(),(),(),(),(22hOyxfyxfyxfphyxfkiiyiiiixii)()()(2hOxyphxyii 將以上結果代入（將以上結果代入（9.149.14）得：）得： )()()(22111kkhxyxyii)()()()()(221hOxyphxyxyhxyiiii )()()()()(32221hOxyphxyhxyiii （9.16） 對式對式(9.15)(9.15)和和(9.16)(9.16)進行比較系數后可知進行比較系數后可知, ,只要只要 211221p（9.17） 成立成立, ,格式格式(9.1

32、4)(9.14)的局部截斷誤差就等于的局部截斷誤差就等于)(3hO有有2 2階階精度精度式式(9.17)(9.17)中具有三個未知量中具有三個未知量, ,但只有兩個方程但只有兩個方程, ,因而因而有無窮多解。若取有無窮多解。若取 , ,則則p p=1=1，這是無窮多解，這是無窮多解中的一個解，將以上所解的值代入式中的一個解，將以上所解的值代入式(9.14)(9.14)并改寫并改寫可得可得 2121),(),()(21121211hkyxfkyxfkkkhyyiiiiii 不難發(fā)現，上面的格式就是改進的歐拉格式。不難發(fā)現，上面的格式就是改進的歐拉格式。凡滿足條件式（凡滿足條件式（9.179.17

33、）有一簇形如上式的計算格式，）有一簇形如上式的計算格式，這些格式統(tǒng)稱為二階龍格這些格式統(tǒng)稱為二階龍格庫塔格式。因此改進的庫塔格式。因此改進的歐拉格式是眾多的二階龍格歐拉格式是眾多的二階龍格庫塔法中的一種特殊庫塔法中的一種特殊格式。格式。 若取若取 , ,則則 ，此時二階龍格，此時二階龍格- -庫塔庫塔法的計算公式為法的計算公式為 0121, 12p)2,(),(1212121khyxfkyxfkhkyyiiiiii1,2 , 1 , 0ni 此計算公式稱為變形的二階龍格此計算公式稱為變形的二階龍格庫塔法。式中庫塔法。式中 為區(qū)間為區(qū)間 的中點。的中點。 21ix1,iixx9.3.3 三階龍格

34、三階龍格- -庫塔法庫塔法 為了進一步提高精度，設除為了進一步提高精度，設除 外再增加一點外再增加一點 pix) 1(qpqhxxiqi并用三個點并用三個點 , , , , 的斜率的斜率k1 1, ,k2 2, ,k3 3加權平均加權平均得出平均斜率得出平均斜率k* *的近似值，這時計算格式具有形式的近似值，這時計算格式具有形式: : ixpixqix11 12233121(,)(,)iiiiiiyyhkkkkf x ykf xph yphk（9.18） 為了預報點為了預報點 的斜率值的斜率值k3 3, ,在區(qū)間在區(qū)間 內有兩內有兩個斜率值個斜率值k1 1和和k2 2可以用可以用, ,可將可將

35、k1 1, ,k2 2加權平均得出加權平均得出 上的平均斜率上的平均斜率, ,從而得到從而得到 的預報值的預報值 qixqiixx,qiixx,)(qixyqiy1122i qiyyqhkk于是可得于是可得 ),(3qiqiyxfk運用運用Taylor展開方法選擇參數展開方法選擇參數 , ,可以使可以使格式格式( (9.18)的局部截斷誤差為的局部截斷誤差為 , ,即具有三階精度，即具有三階精度，這類格式統(tǒng)稱為這類格式統(tǒng)稱為三階龍格三階龍格庫塔方法庫塔方法。下列是其中。下列是其中的一種，稱為的一種，稱為庫塔（庫塔（Kutta）公式。）公式。 12212, ,p q )(4hO)4(6)2(,(

36、)2,(),(3211211312121kkkhyykkhyxfkkhyxfkyxfkiiiiiiii（9.19） 121231232221233132611pqpqpq 1141123122666,1,1,2pq 9.3.4 四階龍格四階龍格庫塔法庫塔法 如果需要再提高精度，用類似上述的處理方法，如果需要再提高精度，用類似上述的處理方法，只需在區(qū)間只需在區(qū)間 上用四個點處的斜率加權平均作為上用四個點處的斜率加權平均作為平均斜率平均斜率k*的近似值，構成一系列四階龍格的近似值，構成一系列四階龍格庫塔公庫塔公式。具有四階精度，即局部截斷誤差是式。具有四階精度，即局部截斷誤差是 。 由于推導復雜，

37、這里從略，只介紹最常用的一種由于推導復雜，這里從略，只介紹最常用的一種四階經典龍格四階經典龍格庫塔公式庫塔公式。 qiixx,)(5hO)22(6),()2,()2,(),(43211314221312121kkkkhyyhkyxfkkhyxfkkhyxfkyxfkiiiiiiiiii（9.20） 9.3.5 9.3.5 四階龍格四階龍格庫塔法算法實現庫塔法算法實現(1)(1) 計算步驟計算步驟 輸入輸入 , ,h,Nh,N 使用龍格使用龍格庫塔公式（庫塔公式（9.209.20）計算出）計算出y y1 1 輸出輸出 ，并使，并使 轉到轉到 直至直至n n N N 結束。結束。 10,xx11,

38、 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 k3) k4 y0+ h (k1+ 2 k2+ + 2 k3+ k4)/6 y1 輸 出 x1, y1 n + 1 n n = N ? x1 x0 y1 y0 結 束 n y (3)(3) 程序實現程序實現( (四階龍格四階龍格- -庫塔法計庫塔法計 算常微分方程初值問題算常

39、微分方程初值問題) ) 例例9.4 9.4 取步長取步長h=0.2h=0.2，用經典格式求解初值問題，用經典格式求解初值問題 1)0(2yxyy10 x解解: : 由四階龍格由四階龍格- -庫塔公式可得庫塔公式可得 2 . 0, 1, 0,2),(00hyxxyyxf0),(001yxfk2 . 0) 1 , 1 . 0()2,(10202fkhyxfkh204. 0)02. 1 , 1 . 0()2,(20203fkhyxfkh41632. 0)0408. 1 , 2 . 0(),(3004fhkyxfkh040811. 1)22(62 . 0432101kkkkyy可同樣進行其余可同樣進行

40、其余y yi i的計算。本例方程的解為的計算。本例方程的解為，從表中看到所求的數值解具有，從表中看到所求的數值解具有4位有效數字。位有效數字。 2xey 龍格龍格庫塔方法的庫塔方法的推導基于推導基于Taylor展開方法，展開方法，因而它要求所求的解具有較好的光滑性。如果解的因而它要求所求的解具有較好的光滑性。如果解的光滑性差，那么，使用四階龍格光滑性差，那么，使用四階龍格庫塔方法求得的庫塔方法求得的數值解，其精度可能反而不如改進的歐拉方法數值解，其精度可能反而不如改進的歐拉方法。在。在實際計算時，應當針對問題的具體特點選擇合適的實際計算時，應當針對問題的具體特點選擇合適的算法。算法。 9.3.

41、6 變步長的龍格變步長的龍格-庫塔法庫塔法 在微分方程的數值解中，選擇適當的步長是非常在微分方程的數值解中，選擇適當的步長是非常重要的。單從每一步看，重要的。單從每一步看，步長越小，截斷誤差就越步長越小，截斷誤差就越小；但隨著步長的縮小，在一定的求解區(qū)間內所要??；但隨著步長的縮小，在一定的求解區(qū)間內所要完成的步數就增加了。這樣會引起計算量的增大，完成的步數就增加了。這樣會引起計算量的增大，并且會引起舍入誤差的大量積累與傳播。因此微分并且會引起舍入誤差的大量積累與傳播。因此微分方程數值解法也有選擇步長的問題。方程數值解法也有選擇步長的問題。 以經典的四階龍格以經典的四階龍格- -庫塔法庫塔法(

42、(9.20)為例。從節(jié)點為例。從節(jié)點x xi i出發(fā)，先以出發(fā)，先以h為步長求出一個近似值，記為為步長求出一個近似值，記為 ，由于局部截斷誤差為由于局部截斷誤差為 ，故有，故有 )(1hiy)(5hO5)(11)(chyxyhii當h h值不大時，式中的系數值不大時，式中的系數c c可近似地看作為常數?？山频乜醋鳛槌?。然后將步長折半然后將步長折半, ,即以為即以為 步長步長, ,從節(jié)點從節(jié)點x xi i出發(fā)出發(fā), ,跨兩跨兩步到節(jié)點步到節(jié)點x xi+1i+1, ,再求得一個近似值再求得一個近似值 , ,每跨一步的截每跨一步的截斷誤差是斷誤差是 , ,因此有因此有2h)2(1hiy52 hc

43、5)2(1122)()(hcxyxyhii這樣這樣 161)()()(11)2(11hiihiiyxyyxy)(151)()(1)2(1)2(11hihihiiyyyxy由此可得由此可得 這表明以這表明以 作為作為 的近似值，其誤差可用步的近似值，其誤差可用步長折半前后兩次計算結果的偏差長折半前后兩次計算結果的偏差 )2(1hiy)(1ixy)(1)2(1hihiyy來判斷所選步長是否適當來判斷所選步長是否適當當要求的數值精度為當要求的數值精度為時：時： （1 1）如果）如果，反復將步長折半進行計算，直，反復將步長折半進行計算，直至至為止為止, ,并取其最后一次步長的計算結果作為并取其最后一次

44、步長的計算結果作為 （2 2）如果）如果為止，并以上一次步長的計算結果作為為止，并以上一次步長的計算結果作為 。 這種通過步長加倍或折半來處理步長的方法稱為這種通過步長加倍或折半來處理步長的方法稱為變步長法。表面上看，為了選擇步長，每一步都要變步長法。表面上看，為了選擇步長，每一步都要反復判斷反復判斷，增加了計算工作量，但在方程的解，增加了計算工作量，但在方程的解y(x)y(x)變化劇烈的情況下，總的計算工作量得到減少，結變化劇烈的情況下，總的計算工作量得到減少，結果還是合算的。果還是合算的。1iy1iy 9.4 算法的穩(wěn)定性及收斂性算法的穩(wěn)定性及收斂性9.4.1 穩(wěn)定性穩(wěn)定性 穩(wěn)定性在微分方

45、程的數值解法中是一個非常穩(wěn)定性在微分方程的數值解法中是一個非常重要的問題。因為微分方程初值問題的數值方法是重要的問題。因為微分方程初值問題的數值方法是用差分格式進行計算的，而在差分方程的求解過程用差分格式進行計算的，而在差分方程的求解過程中，存在著各種計算誤差，這些計算誤差如舍入誤中，存在著各種計算誤差，這些計算誤差如舍入誤差等引起的擾動，在傳播過程中，可能會大量積累，差等引起的擾動，在傳播過程中，可能會大量積累，對計算結果的準確性將產生影響。這就涉及到算法對計算結果的準確性將產生影響。這就涉及到算法穩(wěn)定性問題。穩(wěn)定性問題。 當在某節(jié)點上當在某節(jié)點上x xi i的的y yi i值有大小為值有大

46、小為的擾動時，的擾動時，如果在其后的各節(jié)點如果在其后的各節(jié)點 上的值上的值y yi i產生的偏差都產生的偏差都不大于不大于，則稱這種方法是穩(wěn)定的。，則稱這種方法是穩(wěn)定的。 穩(wěn)定性不僅與算法有關，而且與方程中函數穩(wěn)定性不僅與算法有關，而且與方程中函數f(x,y)f(x,y)也有關，討論起來比較復雜。也有關，討論起來比較復雜。為簡單起見，為簡單起見，通常只針對模型方程通常只針對模型方程)(ijxj)0(yy來討論。一般方程若局部線性化，也可化為上述形來討論。一般方程若局部線性化，也可化為上述形式。模型方程相對比較簡單，若一個數值方法對模式。模型方程相對比較簡單，若一個數值方法對模型方程是穩(wěn)定的，并

47、不能保證該方法對任何方程都型方程是穩(wěn)定的，并不能保證該方法對任何方程都穩(wěn)定，但若某方法對模型方程都不穩(wěn)定，也就很難穩(wěn)定，但若某方法對模型方程都不穩(wěn)定，也就很難用于其他方程的求解。用于其他方程的求解。先考察顯式先考察顯式EulerEuler方法的穩(wěn)定性。模型方程方法的穩(wěn)定性。模型方程的的EulerEuler公式為公式為 yy)(),(1iiiiiiyhyyxhfyy), 2 , 1 , 0()1 (iyhi將上式反復遞推后，可得將上式反復遞推后，可得 011)1(yhyii), 2 , 1()1 (00iyyhyiii或或h1式中式中 要使要使y yi i有界，其充要條件為有界，其充要條件為 1

48、11h即即 由于由于 ，故有，故有 020 h（9.269.26） 可見，如欲保證算法的穩(wěn)定，顯式可見，如欲保證算法的穩(wěn)定，顯式EulerEuler格式格式的步長的步長h h的選取要受到式（的選取要受到式（9.269.26）的限制。）的限制。 的絕的絕對值越大，則限制的對值越大，則限制的h h值就越小。值就越小。 用隱式用隱式EulerEuler格式，對模型方程格式，對模型方程 的計算公式為，可化為的計算公式為，可化為)(11iiiyhyyiiyhy111由于由于 , ,則恒有則恒有 , ,故恒有故恒有 0111hiiyy1 因此，隱式因此，隱式EulerEuler格式是絕對穩(wěn)定的（無條件穩(wěn)格

49、式是絕對穩(wěn)定的（無條件穩(wěn)定）的（對任何定）的（對任何h0h0）。）。 9.4.2 9.4.2 收斂性收斂性 常微分方程初值問題的求解，是將微分方程轉化常微分方程初值問題的求解，是將微分方程轉化為差分方程來求解，并用計算值為差分方程來求解，并用計算值y yi i來近似替代來近似替代y(xy(xi i) )，這種近似替代是否合理，還須看分割區(qū)間這種近似替代是否合理，還須看分割區(qū)間 的長度的長度h h越來越越來越小時，即小時，即 時，時， 是否成立。若是否成立。若成立，則稱該方法是收斂的，否則稱為不收斂。成立，則稱該方法是收斂的，否則稱為不收斂。 這里仍以這里仍以EulerEuler方法為例，來分析

50、其收斂性。方法為例，來分析其收斂性。EulerEuler格式格式iixx,101iixxh)(iixyy),(1iiiiyxhfyy設設 表示取表示取 時時, 按按Euler公式的計算結果公式的計算結果, 即即1iy)(iixyy )(,()(1iiiixyxhfxyyEuler方法局部截斷誤差為：方法局部截斷誤差為：)()(2)(1211 iiiixxyhyxy設有常數設有常數 ，則，則 )(max21xycbxa 211)(chyxyii（9.27） 總體截斷誤差總體截斷誤差 1111111)()(iiiiiiiyyyxyyxy)(,()(),(11iiiiiiiixyxhfxyyxhfy

51、yy),()(,()(iiiiiiyxfxyxfhyxy（9.28） 又又 由于由于f(x,y)f(x,y)關于關于y y滿足李普希茨條件滿足李普希茨條件, ,即即 iiiiiiyxyLyxfxyxf)(),()(代入代入(9.28)上式，有上式，有 iiiiyxyhLyy)()1(11ihL)1( 再利用式（再利用式（9.27），式（），式（9.28） 1111111)()(iiiiiiiyyyxyyxy2)1 (chhLi21)1 (chhLii即即 上式反復遞推后，可得上式反復遞推后，可得 1020)1 ()1 (ikkiihLchhL1)1 ()1 (0iihLLchhL（9.29）

52、設設 （T T為常數）為常數） Tihxxi0hLehL 1TLihLieehL)1 (因為因為 所以所以 把上式代入式（把上式代入式（9.299.29），得），得 ) 1(0TLTLieLche若不計初值誤差，即若不計初值誤差，即 ，則有，則有 00heLcTLi)1(（9.30） 式式(9.30)(9.30)說明說明, ,當時當時 , , , ,即即 , ,所以所以EulerEuler方法是收斂的，且其收斂速度為方法是收斂的，且其收斂速度為 ，即具，即具有一階收斂速度。同時還說明有一階收斂速度。同時還說明EulerEuler方法的整體截斷方法的整體截斷誤差為誤差為 ，因此算法的精度為一階。

53、，因此算法的精度為一階。 0h0i)(iixyy )(hO)(hO9.5 亞當姆斯方法亞當姆斯方法9.5.1 亞當姆斯格式亞當姆斯格式 龍格龍格-庫塔方法是一類重要算法，但這類庫塔方法是一類重要算法，但這類算法在每一步都需要先預報幾個點上的斜率算法在每一步都需要先預報幾個點上的斜率值，計算量比較大?？紤]到計算值，計算量比較大?？紤]到計算yi+1之前已得之前已得出一系列節(jié)點上出一系列節(jié)點上 的斜率值，的斜率值，能否利用這些已知值來減少計算量呢？能否利用這些已知值來減少計算量呢？ 這就是亞當姆斯（這就是亞當姆斯（Adams）方法的設計）方法的設計思想。思想。 11,xxxii 設用設用x xi i

54、,x,xi+1i+1兩點的斜率值加權平均作為區(qū)間兩點的斜率值加權平均作為區(qū)間 上的平均斜率，有計算格式上的平均斜率，有計算格式 1,iixx),(),()1 (11111iiiiiiiiiiyxfyyxfyyyhyy（9.21） 選取參數選取參數，使格式（，使格式（9.219.21）具有二階精度。）具有二階精度。 將將 在在x xi i處處Taylor展開展開 1iy)()(! 21)(321hOhyhyyyiii 代入計算格式代入計算格式(9.21)(9.21)化簡化簡, ,并假設并假設, , 因此有因此有 )(),(11iiiixyyxyy)()()()(321hOxyhxyhxyyiii

55、i 與與y(xi+1)在在xi處的處的Taylor展開式展開式 )()(21)()()(321hOxyhxyhxyxyiiii 相比較相比較, 需取需取 21才使格式才使格式(9.21)具有二階精度。這樣導出的計算格式具有二階精度。這樣導出的計算格式 )3(211iiiiyyhyy稱之為二階亞當姆斯格式。類似地可以導出三階亞稱之為二階亞當姆斯格式。類似地可以導出三階亞當姆斯格式當姆斯格式。 )51623(12211iiiiiyyyhyy和和四階亞當姆斯格式。四階亞當姆斯格式。 )9375955(243211iiiiiiyyyyhyy( (9.22) 這里和下面均記這里和下面均記 ),(kiki

56、kiyxfy 上述幾種亞當姆斯格式都是顯式的，算法比較上述幾種亞當姆斯格式都是顯式的，算法比較簡單，但用節(jié)點簡單，但用節(jié)點 的斜率值來預報的斜率值來預報區(qū)間區(qū)間 上的平均斜率是個上的平均斜率是個外推過程，外推過程，效果不夠效果不夠理想。為了進一步改善精度，變外推為理想。為了進一步改善精度，變外推為內插內插，即，即增加節(jié)點增加節(jié)點xi+1的斜率值來得出的斜率值來得出 上的平均斜率。上的平均斜率。譬如考察形如譬如考察形如 11,xxxii1,iixx1,iixx),(),()1 (11111iiiiiiiiiiyxfyyxfyyyhyy（9.23） 的隱式格式的隱式格式, ,設設( (9.23)右

57、端的右端的Taylor展開有展開有 )(),(11iiiixyyxyy)()()1 ()()(321hOxyhxyhxyyiiii 可見要使格式可見要使格式(9.23)(9.23)具有二階精度具有二階精度, ,需令需令 , ,這樣就可構造二階隱式亞當姆斯格式這樣就可構造二階隱式亞當姆斯格式 21)(211iiiiyyhyy其實是梯形格式。類似可導出三階隱式亞當姆斯格其實是梯形格式。類似可導出三階隱式亞當姆斯格式式 )85(12111iiiiiyyyhyy和四階隱式亞當姆斯格式和四階隱式亞當姆斯格式 )5199(242111iiiiiiyyyyhyy（9.249.24） 9.5.2 亞當姆斯預報

58、亞當姆斯預報-校正格式校正格式 參照改進的歐拉格式的構造方法，以四階亞當參照改進的歐拉格式的構造方法，以四階亞當姆斯為例，將顯式（姆斯為例，將顯式（9.22）和隱式（）和隱式（9.24）相結合，）相結合，用顯式公式做預報，再用隱式公式做校正，可構成用顯式公式做預報，再用隱式公式做校正，可構成亞當姆斯預報亞當姆斯預報-校正格式校正格式 )9375955(243211iiiiiiyyyyhyy),(111iiiyxfy)5199 (242111iiiiiiyyyyhyy),(111iiiyxfy（9.259.25） 預報：預報： 校正：校正： 這種預報這種預報- -校正格式是四步法，它在計校正格式

59、是四步法，它在計算算y yi+1i+1時不但用到前一步的信息時不但用到前一步的信息 ，而且，而且要用到再前面三步的信息要用到再前面三步的信息 ，因，因此它不能自行啟動。在實際計算時，可借助此它不能自行啟動。在實際計算時，可借助于某種單步法，譬如四階龍格于某種單步法，譬如四階龍格庫塔法提供庫塔法提供開始值開始值 。 iiyy,321,iiiyyy321,yyy例例9.5 取步長取步長h=0.1,h=0.1,用亞當姆斯預報用亞當姆斯預報- -校正公式求解校正公式求解 初值問題初值問題 1)0(yyxy10 x的數值解。的數值解。 解解: : 用四階龍格用四階龍格- -庫塔公式求出發(fā)值庫塔公式求出發(fā)

60、值 , ,計算得：計算得：1 . 0, 1, 0,),(00hyxyxyxf321,yyy399717. 1,242805. 1,110342. 1321yyy表中的表中的 ,y,yi i和和y(xy(xi i) )分別為預報值、校正值和準確分別為預報值、校正值和準確解解( ),( ),以比較計算結果的精度。以比較計算結果的精度。 再使用亞當姆斯預報再使用亞當姆斯預報- -校正公式校正公式(9.25),(9.25),iy12xeyx9.6 一階常微分方程組與高階方程一階常微分方程組與高階方程 我們已介紹了一階常微分方程初值問題的各種我們已介紹了一階常微分方程初值問題的各種數值解法，對于一階常微