1、第五章 常微分方程數(shù)值解法1 引言在實際工作中常常會遇到求解常微分方程的定解問題。雖然我們已經(jīng)學(xué)過一些常微分方程的定解問題的解法，但是那只是對一些特殊類型的方程給出了解析方法。而實際問題中歸納出來的常微分方程往往不能用解析方法來求解。有時雖然能求出其解析式子，但真要求它在某一點的值時，還會遇到一些麻煩。見下例。y 1 2xy例 已知初值問題求其解。y(0) 02x 2易知， y(x) e x 0 e dx。若要求它在某一點的值，還要求積分，而這個積分是不能求出的，需要用數(shù)值積分來求解。因此不如直接用數(shù)值解法求其解了。以下我們就研究常微分方程的數(shù)值解法。本章主要考慮下列一階方程的定解問題(又稱初

2、值問題)。1-1 )y f ( x, y)y(x0) y0這里假設(shè)函數(shù)f(x, y)滿足Lipschitz(李卜西茲)條件。 即 f (x, y)滿足：f (x, y) f (x, y) L y y其中， L 為某一常數(shù)。這是保證問題(1-1)有唯一解。下面就研究關(guān)于問題(1-1)的數(shù)值解法。所謂數(shù)值解法，就是直接尋求解函數(shù)y(x)在一系列離散點x1x2xnxn 1上的近似值y1，y2，yn，yn1，。其中，xi1xihi稱為步長，今后如果不特殊說明，我們總是假定是等步長的。即有xi 1 xi h ， 此 時 節(jié) 點 xn x0 nh ， n 0,1,2, 。 由 于 y(x0) y0為已知，

3、所以自然設(shè)想利用這個已知信息求出y(x1)的近似值y1 ，然后由y1求得y(x2)的近似值y2，如此繼續(xù)下去，這就是初值問題數(shù)值解法的一般思想。稱為“步進法”。§ 2 Euler 方法(折線法)2-1 Euler 公式這一方法是初值問題數(shù)值解中最簡單的一個方法，其精度不高。 所以實際計算中很少應(yīng)用。但在某種程度上反映了數(shù)值方法的基本思想，且在此基礎(chǔ)上得到的某些改進的方法目前還在使用，因此我們有必要介紹一下這種方法。Euler 公式的推導(dǎo)方法很多，在此我們用Tarlor 展開的方法。h2y(xn 1)y(xnh)y(xn)y (xn)hy ( n)2!其中， n(xn,xn 1) ，由

4、(1-1)式知，h2y(xn 1) y(xn) hfxn, y(xn)2! y ( n)n 0,1,2,h2當 h 充分小時，略去h y ( n ) ，即得2! ny(xn 1) y(xn) hfxn, y(xn)取近似，并寫成等式得，2-1 )yn 1yn hf(xn, yn) n 0,1,2,此即稱為Euler 公式。h2由于略去的是h y ( n)，可見這就是誤差。但這只是在計2! n算第 n 步時產(chǎn)生的誤差，并非是從一開始到現(xiàn)在所產(chǎn)生的誤差。因此我們稱這個誤差為局部截斷誤差。即在假定yny(xn )的假設(shè)下，計算y(xn 1)時產(chǎn)生的誤差y(xn 1) yn 1，稱為局部截斷誤h2差。

5、由此知Euler 公式的局部截斷誤差約為h y (xn )。2! n2-2 后退的 Euler 公式后退的 Euler 公式的推導(dǎo)同樣也有許多方法。在此，我們用差商代替導(dǎo)數(shù)的方法。因為y(xn1)y(xn 1) y(xn)xn 1 xn假設(shè)yny(xn)，又有(1-1)知y(xn 1)f(xn 1,y(xn 1)，取近似整理得，yn 1 yn hf(xn 1,yn 1) n 0,1,2,( 2-2)此即稱為后退的(隱式)Euler 公式。公式( 2-2)在計算yn 1時要用yn 1。這樣的公式稱為隱式公式，而象公式(2-1 )的式子稱為顯式公式。對于隱式公式，在計算時要先給yn 1 提供一個初

6、值，然后再用給定的公式開始計算，通常稱為迭代法。對(2-2)式的迭代公式如下，k 0,1,2-3)yn(01)ynhf(xn,yn)yn(k11)ynhf(xn 1,yn(k1)直到y(tǒng)n(k11)yn(k1)為止。以下研究式(2-2)的局部截斷誤差。 yn 1ynhf(xn 1, yn 1) yn hy(xn 1)又 y (xn 1 ) y (xn h) y (xn) y ( n )h 代入上式得，yn 1yn hy(xn) y ( n)h2h2又有y(xn 1) y(xn) hf(xn, y(xn)2! y ( n) 下式減上式，且注意到已假設(shè)yn y(xn)，故得，y(xn 1) yn 1

7、= h y ( n)2h2后退的Euler 公式的局部截斷誤差約為h y (xn)。2! n2-3 梯形公式今考察以上兩個公式，可見它們的局部截斷誤差相差一個負號，即有h23yn 1yn hf(xn,yn) 2! y (xn) o(h )h23yn 1yn hf(xn 1, yn 1) y (xn) o(h )2!h2以上兩式相加除以2，則可以消去h y (xn )，從而得到2! nh3yn 1yn 2 f(xn,yn) f(xn 1, yn 1) o(h )略去o(h3)，則得一公式，hyn 1yn 2f(xn,yn) f(xn 1,yn 1)( 2-4) 此稱為梯形公式。其幾何意義如圖。注

8、意到(2-4)式為隱式公式，同樣需要迭代法來求解。其迭代公式為，(0)yn 1ynhf(xn,yn)(k 1) h(k) k 0,1,( 2-5)yn 1yn2f (xn,yn)f (xn 1,yn 1)以下分析迭代公式(2-5 )的收斂性。( 2-4)減去(2-5)得yn 1 yn(k11)h f(xn 1,yn 1) f (xn 1,yn(k)1)2由于已假設(shè)f(x, y)滿足Lipschitz(李卜西茲)條件。 即 f(x, y)滿足：f (x, y) f (x, y) L y y所以有yn1 yn(k11)L2h yn1yn(k1)，取L2h1，則有yn 1yn(k11)yn 1 yn

9、( k1)從而收斂(帶條件收斂)。2-4 改進的 Euler 公式由(2-5)式知，用梯形公式計算時要反復(fù)求函數(shù)值，因而工作量是很大的。而改進的Euler 公式是指在用梯形公式迭代時只迭代一次便取作yn 1 ，這就是，預(yù)測：yn 1ynhf (xn, yn )h校正：yn1yn2 f(xn,yn)f(xn1,yn1)( 2-6)以上稱為預(yù)測-校正系統(tǒng)。也稱為改進的Euler公式為便于上機常采用以下形式，2-7)yp yn hf(xn,yn)yc yn hf (xn 1, yp)1yn 1(yp yc)例 (見教材)在以上推導(dǎo)后退的Euler 公式時，我們用差商代替了導(dǎo)數(shù)。但若改用中心差商代替導(dǎo)

10、數(shù)就有y(xn1)2hy(xn1) y(xn)取近似并整理得，y(xn 1) y(xn 1) 2hf (xn, yn)從而得，yn 1 yn 1 2hf (xn , yn )用此公式作為預(yù)測值，用梯形公式作為校正值，便得到以下的預(yù)測 -校正系統(tǒng)，預(yù)測：yn 1 yn 1 2hf(xn,yn)h校正：yn 1 yn f(xn,yn) f(xn 1, yn 1)( 2-8)( 2-8)與(2-6)相比，突出的特點是它們有相同的精度。不過( 2-8)式是兩步公式，而前面幾式是單步公式。兩步以上的公式稱為多步法，多步法不能自開始(自起步)，需要借助其他單步法才能開始進行。而單步法可以自開始。§

11、; 3 Runge-Kutta 方法在§ 1 中， 我們曾經(jīng)用Taylor 展開的方法推導(dǎo)了Euler 公式。那時我們略去了o(h2)項，自然想到如果多取幾項，也就是略去o(hp 1)項( p取的更大些)應(yīng)該得到精度更高的公式。這就是Taylor 級數(shù)法。它的一般公式如下，h2hpyn 1ynhynynyn(p)( 3-1)( 3-1)式的局部截斷誤差為hp1y(xn 1)yn 1=y(p 1)( n) n(xn,xn 1)(p 1)!一般地有以下定義，定義 若一種方法的局部截斷誤差為o(hp 1)，則稱該方法具有 p 階精度。由( 3-1)知，當p=1 時， ( 3-1 )式即為E

12、uler 公式。顯然Euler公式具有1 階精度。同理，后退的Euler 公式也具有1 階精度；而梯形公式具有2 階精度。由( 3-1)知， p 越大精度越高。但是它要用到f(x,y)的高階導(dǎo)數(shù)， 而求 f(x,y)的高階導(dǎo)數(shù)是不容易求出的。因此， Tarlor 級數(shù)法在實際問題中很少應(yīng)用，它的作用在于啟發(fā)我們?nèi)ヌ剿鞲玫姆椒ā?-1 Runge-Kutta 方法的基本思想根據(jù)微分中值定理，存在01 ，使得y(xn 1) y(xn)y(x h)hn由方程(1-1)知，y (xnh) f (xn h, y(xnh)有y(xn 1) y(xn) h f(xnh,y(xnh)( 3-2)令K f (

13、xnh,y(xnh) 稱為區(qū)間xn,xn 1 上的平均斜率。由此可見，只要對這個平均斜率提供一種算法，由(3-2)式便可以得到一種計算公式。若取0，即取(xn, yn)點的斜率作為整個區(qū)間xn,xn 1上的平均斜率，這時K* f (xn, yn)，則得Euler 公式y(tǒng)n 1yn hf (xn , yn) n 0,1,2,若取1，即取(xn 1, yn 1)點的斜率作為整個區(qū)間xn,xn 1上的平均斜率，這時K* f (xn 1, yn 1)，則得后退的Euler 公式y(tǒng)n1yn hf(xn 1,yn 1) n 0,1,2,若取K1f(xn,yn)，K2f(xn1,ynhf (xn,yn)，令

14、*1K*1 (K1 K2)2即取(xn, yn) 和 (xn 1, yn 1) 兩點斜率的算術(shù)平均作為整個區(qū)間xn,xn 1上的平均斜率，則得下列公式，yn 1yn h (K1 K2)2K1f(xn,yn)( 3-3)K2f(xn1,yn hK1)這顯然是改進的Euler 公式(即梯形公式)。我們知道，改進的Euler 公式是具有2 階精度，而Euler 公式和后退的Euler 公式只有1 階精度?？梢娪脙蓚€點斜率的算術(shù)平均作為平均斜率比只取一個點的斜率作為平均斜率精度要高些。由此得，如果設(shè)法在xn,xn 1內(nèi)多預(yù)測幾個點的斜率值，然后將它們加權(quán)平均作為平均斜率K *，則有可能構(gòu)造出具有更高精

15、度的計算公式。這就是Runge-Kutta 方法的基本思想。3-2 二階 Runge-Kutta 方法公式 ( 3-3) 是一個特殊的二階Runge-Kutta 方法。 它用了xn，xn 1這兩個點上的斜率值。我們現(xiàn)在在xn,xn 1內(nèi)取xn， xn p 0 p 1將這兩個點上的斜率值K1， K 2的線性組合作為平均斜率K*。 顯然， K1 f (xn,yn)， 而 K2f (xn p, yn p)， 對于其中的yn p我們用 Euler 公式來計算得，yn pynphf (xn , yn)( 3-4)為此有，yn 1 yn h( 1K12K2)K1f(xn,yn)( 3-5)K2f (xn

16、p, ynphK1)其中含有三個待定參數(shù)1, 2, p。確定這三個待定參數(shù)的原則是使得公式 ( 3-5) 具有二階精度。為此將 K2在 (xn, yn )點作二元函數(shù) Taylor 展開，得K2f(xn,yn) ph(fx f fy)n將 K1， K2代入(3-5)的第一式得yn 1yn h( 1k12K2)yn h 1fn 2fn2ph(fx f fy)nyn ( 12)hfn2ph2(fxf fy)n而二階 Taylor 級數(shù)法公式為h2yn 1 yn hynyn1, 2, p應(yīng)滿足3-6)12112p 2滿足（3-6）形如（3-5）的所有公式均稱為二階Runge-Kutta 公1式。特別

17、當121 ， p=1 時就是改進的Euler 公式。3-3 三階 Runge-Kutta 方法為提高精度，我們現(xiàn)在在xn,xn 1內(nèi)取xn，xnp， xn q0p q 1三個點，用這三個點上的斜率值K1 ， K2， K 3的線性組合作為平均斜率 K*。顯然，K1 f(xn,yn)，K2 f(xn p,ynphK1)，而K 3 的預(yù)測值我們用K3 f(xn q,yn q)f (xn qh, yn qh(rK1 sK2)來計算。從而得，yn 1 yn h( 1K12 K23K3)K1f(xn,yn)( 3-7)K2 f(xn ph, yn phK1)K3 f(xn qh, yn qh(rK1 sK

18、2)為了確定（3-7）中的待定參數(shù)1, 2, 3, p,q,r,s。我們同樣用將K1， K2， K3作 Taylor 展開，然后代入（3-7）中的yn 1。通過與 三階 Taylor 級數(shù)法比較系數(shù)，則得待定參數(shù)滿足的條件為r s 1, 1231( 3-8)12p 3q 21以及2 p23q2( 3-9)23133pqs36參數(shù)滿足 ( 3-8)( 3-9) 的形如 ( 3-7) 的公式統(tǒng)稱為三階Runge-Kutta公式。它也是一族公式。用同樣的方法我們還可以得到四階Runge-Kutta 公式。Runge-Kutta 公式的一般表示為ryn 1 yn h iKii1K1f(xn,yn)i1

19、Kif(xnpih, yn hrijKj)i 2,3, ,rj1取r=4,用同樣的方法即可以得到四階Runge-Kutta 公式。下面就是一個經(jīng)典的四階Runge-Kutta 公式。yn 1 yn h(K1 2K2 2K3 K4)6K1f (xn, yn )K2f(xn h,ynhK1)( 3-10)22K3f(xn2, yn 2 K2)K4 f(xn h,yn hK3)可以證明其截斷誤差為o(h5 ) 。應(yīng)注意的是，Runge-Kutta 公式的推導(dǎo)基于Taylor展開方法，所以它要求解函數(shù)有很好的光滑性。即y(x) 要具有所要求的導(dǎo)數(shù) 。 若 不 然 ， 用 高 階 Runge-Kutta

20、 公 式 可 能 不 如 用 低 階 Runge-Kutta 公式效果好。§ 4 單步法的收斂性與穩(wěn)定性4-1 收斂性定義 2 若一種數(shù)值方法對于任意固定的xnx0 nh，單步法yn 1 yn h (xn, yn,h)( 4-1)滿足ynh 0(n ) y(xn)，則稱該方法是收斂的。yy例 考慮( 4-2)y(0)y0解 ( 4-2)的Euler 公式為yn 1yn h yn (1h)yn，所以有yn (1 h)ny0 ，因為 x0 0，xnnh，11yn(1 h) h hny0 (1 h) h xny0 lim yny0e xny(xn)。h0對于一般單步法yn 1ynh (xn

21、, yn,h)( 4-1)有以下結(jié)論，定理 假設(shè)單步法有p 階精度， 且增量函數(shù)(x, y,h)關(guān)于y 滿足 Lipschitz 條件。即對任意y, y，總存在常數(shù)L>0, 使(x, y,h) (x, y,h) L y y( 4-3)則單步法(4-1 )收斂。且整體截斷誤差為y(xn) yn o(hp)。定理的證明從略。作為應(yīng)用我們考慮梯形公式(改進的Euler方法)的收斂性。梯形公式的增量函數(shù)為(x, y,h)12 f (x, y)f(xh,yhf (x, y)而1(x, y,h)2 f (x, y)f(xh,yhf (x, y)以上兩式相減，并注意到不等式的性質(zhì)，得1(x,y,h)(

22、x,y,h) 2 f(x,y) f(x,y)12 f(x h,y hf (x,y) f (x h,y hf(x, y)又由于我們已假設(shè)函數(shù)f(x,y) 關(guān)于y 滿足 Lipschitz 條件，所以上式變?yōu)?x, y,h) (x,y,h)3 2 y yL2 y hf (x, y) y hf (x, y)整理得，(x,y,h)(x,y,h)L LhL2y y L y y222所以，梯形公式的增量函數(shù)關(guān)于y 滿足 Lipschitz 條件。從而收斂。4-2 穩(wěn)定性定義 3 若一種數(shù)值方法在節(jié)點xn上的值yn有大小的擾動，而以后節(jié)點xm (m>n)上的值ym所產(chǎn)生的擾動不超過，則稱 該方法是穩(wěn)定

23、的。由于穩(wěn)定性的討論比較復(fù)雜，所以通常對模型方程y y (0 )( 4-4)進行討論，以得到穩(wěn)定性條件。所以穩(wěn)定性的定義還可以如下描述，定義 4 用單步法(4-1)解模型方程(4-4)所得到的穩(wěn)定性方程 n 1 E( h) n，若 E( h) 1，則稱此法是絕對穩(wěn)定的。由 E( h) 1，所得區(qū)間稱為穩(wěn)定區(qū)間。n為在節(jié)點值yn上的擾動值。今考慮用Euler 公式解模型方程(4-4)時的穩(wěn)定區(qū)間。為此設(shè) yn 為理論值，yn* 為實際計算值，則有理論值yn 1yn h yn (1 h)yn實際計算值yn 1(1h)yn兩式相減得穩(wěn)定性方程，n 1 (1 h) n 其中n為第n 次擾動。顯然，要保

24、證穩(wěn)定的條件是n 1 n ，即有 1 h 1 。解之得 Euler 公式的穩(wěn)定區(qū)間為2 h0對于后退的Euler 公式，則有yn 1yn h yn 1*yn 1yn h yn 1兩式相減得穩(wěn)定性方程，n 1 n h n 1 ，解得 n 1n1h1由于0，對任意步長h>0 總有 11，所以我們通常說1h后退的 Euler 公式是無條件穩(wěn)定的?；蛘哒f后退的Euler 公式的穩(wěn)定區(qū)間為h 0 。同樣方法我們可以討論其他公式的穩(wěn)定區(qū)間。比如， 改進的Euler 公式(2-6) 。由于公式為預(yù)測：yn1ynhf (xn,yn)h校正：yn 1yn f (xn,yn)f (xn 1,yn 1)2對于

25、模型(4-4)變?yōu)轭A(yù)測：yn 1ynhyn校正：yn 1yn2h ynyn1h即yn 1 yn yn( ynhyn)2ynh2 ynyn2hyn1 221 h h y2n1yn* 11 h 12 2h2yn*1兩式相減得穩(wěn)定性方程，n 1 1 h 1 2h2 n 得穩(wěn)定性 條件為其他公式不再討論。例 見教材5 線性多步法線性多步法的公式為1 h 1 2h2 122 h0這一節(jié)我們研究線性多步法。解之得穩(wěn)定區(qū)間為5-1 )rryn 1kyn k h kyn kk0k 1其中，yn k f （xn k,yn k）。當 1 0， （ 5-1）為顯式公式；當0， （ 5-1）為隱式公式。適當選擇（5-

26、1）式中的待定參數(shù)，就可以得到一系列的線性多步法公式。確定待定參數(shù)時常用的方法有數(shù)值積分法和Taylor 展開法，在此我們用Taylor 展開法來 確定待定參數(shù)，因為這種方法更具有一般性，且更容易理解。設(shè) yn ky(xn k ) y(xn kh)yn ky (xn k) y (xn kh) 分別在xn點作Taylor 展開，得， yn kyn kpj0pj1( kh)jj!( kh)(j)( kh)p (p 1)yn(p 1)! ynj1(j 1)!(j)( kh)p (p 1)ynynp!5-1）并整理得ryn 1 (k)ynk0pjrrh ( k)j k j ( k)j 1 k yn(j

27、)j 1 j! k 1k 1hp1 rrh ( k)p 1 k (p 1)( k)p k yn(p 1)( p 1)! k 1k 1要使上式具有p 階精度， 即局部截斷誤差為o(hp 1)。 又由y(xn 1)Taylor 展開式p hj y(xn1) j0 j!ynhp1(p 1)!yn(p 1)知，需要符合到h p 項，所以應(yīng)有下式成立r5-2)k1k0rr( k)j k j ( k)j1 k 1k1k 15-2)式成立時，局部截斷誤差為p1rry(xn 1) yn 1 h 1( k)p1 k (p 1)( k)p k yn(p1)(p 1)! k 1k 1( 5-3)以下我們具體的考慮兩

28、個四步方法。先考慮顯式公式y(tǒng)n 10yn1yn 12yn 2h( 0yn 1yn 12yn 23yn 3)要使上式有四階精度，根據(jù)(5-2)式知其系數(shù)應(yīng)滿足0121215-4)120123142214263131 8 2 3 112 2 27 31121231 16 2 4 1 32 2 108 3 17 個待定參數(shù)，而只有5 個方程，所以（5-4）的解是不唯一的。若取120，則得10,055593724, 124, 2 24924此時的公式為hyn 1 yn(55fn 59fn 1 37fn 2 9fn 3)( 5-5)24其中， fn k f(xn k, yn k)。 ( 5-5)稱為四步

29、Adams 顯式公式。再考慮隱式公式（用yn 1 代替yn 3） ，即yn 10yn1yn 12yn 2h( 1yn 10yn1yn 12yn 2)要使上式有四階精度，根據(jù)（5-2）式知其系數(shù)應(yīng)滿足10125-6)1221012114221214211211218 23 13 112 211 16 2 4 1 4 1 32 2 17 個待定參數(shù)，而只有5 個方程，所以（5-6）的解是不120 ，則得01,1919524, 0 24, 124124此時的公式為hyn 1 yn (9fn 1 15fn 5fn1 fn 2)( 5-7)24其中，fn kf(xn k, yn k)。 ( 5-7)稱為

30、四步Adams 隱式公式。用以上的方法還可以得到Milne 公式以及Hamming 公式等。作為以上方法的應(yīng)用我們考慮下例，y f ( x, y)例 解初值問題用顯式二步公式y(tǒng)(x0)y0yn 10yn1yn 1 h( 0yn1yn 1)。( 5-8)試確定待定參數(shù)0, 1，0, 1使方法階數(shù)盡可能高；并求局部截斷誤差。解 由于以上公式不易記憶，這里我們直接用Taylor 展開方法。因為yn 1y(xnh)ynhynhynhynhyn( 4)o(h5)h63 yn(4)o(h4 )h2yn 1 y (xn h) yn hyn yn2其中，yn y(xn) f (xn,yn)代入(5-8)式得，

31、yn 10yn1(ynhynh ynh yn h yn(4) o(h5)23h 0yn1(ynhynh yn h yn(4) o(h4)261( 01)yn (101)hyn (11)h2yn21111( 4 (4)5h yn o(h )6111)h3yn(1111)h4yn(4)o(h5)624!6考慮 Taylor 級數(shù)法公式h2h3hp(p)yn 1ynhynyn ynyn2!3!p!為使方法階數(shù)盡量高，需有011112 1解之得，04 ，1 5，0 4，12。此時公式為三階。所得二步法為yn 14yn 5yn1 2h(2yn yn 1)11將 0, 1，0, 1代入 ( 1 1 1 1)h4yn(4) o(h5)得4!61 4 ( 4)58 h yno(h )1而 Taylor 級數(shù)法公式當p=3 時相應(yīng)的項為h yn o(h )4! n( 11)h4yn(4) o(h5)4!8作為練習(xí)請你用以上方法解下例：y f ( x, y)例 解初值問題用顯式二步公式y(tǒng)(x0 ) y0yn 1(yn yn 1) h( 0yn1yn 1)。試確定待定參數(shù), ，0, 1使方法階數(shù)盡可能高；并求局部截