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

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

3、所以自然設想利用這個已知信息求出y(x1)的近似值y1 ，然后由y1求得y(x2)的近似值y2，如此繼續下去，這就是初值問題數值解法的一般思想。稱為“步進法”。§ 2 Euler 方法(折線法)2-1 Euler 公式這一方法是初值問題數值解中最簡單的一個方法，其精度不高。 所以實際計算中很少應用。但在某種程度上反映了數值方法的基本思想，且在此基礎上得到的某些改進的方法目前還在使用，因此我們有必要介紹一下這種方法。Euler 公式的推導方法很多，在此我們用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 步時產生的誤差，并非是從一開始到現在所產生的誤差。因此我們稱這個誤差為局部截斷誤差。即在假定yny(xn )的假設下，計算y(xn 1)時產生的誤差y(xn 1) yn 1，稱為局部截斷誤h2差。

5、由此知Euler 公式的局部截斷誤差約為h y (xn )。2! n2-2 后退的 Euler 公式后退的 Euler 公式的推導同樣也有許多方法。在此，我們用差商代替導數的方法。因為y(xn1)y(xn 1) y(xn)xn 1 xn假設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)直到yn(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) 下式減上式，且注意到已假設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由于已假設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)式知，用梯形公式計算時要反復求函數值，因而工作量是很大的。而改進的Euler 公式是指在用梯形公式迭代時只迭代一次便取作yn 1 ，這就是，預測：yn 1ynhf (xn, yn )h校正：yn1yn2 f(xn,yn)f(xn1,yn1)( 2-6)以上稱為預測-校正系統。也稱為改進的Euler公式為便于上機常采用以下形式，2-7)yp yn hf(xn,yn)yc yn hf (xn 1, yp)1yn 1(yp yc)例 (見教材)在以上推導后退的Euler 公式時，我們用差商代替了導數。但若改用中心差商代替導

10、數就有y(xn1)2hy(xn1) y(xn)取近似并整理得，y(xn 1) y(xn 1) 2hf (xn, yn)從而得，yn 1 yn 1 2hf (xn , yn )用此公式作為預測值，用梯形公式作為校正值，便得到以下的預測 -校正系統，預測：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 中， 我們曾經用Taylor 展開的方法推導了Euler 公式。那時我們略去了o(h2)項，自然想到如果多取幾項，也就是略去o(hp 1)項( p取的更大些)應該得到精度更高的公式。這就是Taylor 級數法。它的一般公式如下，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)的高階導數， 而求 f(x,y)的高階導數是不容易求出的。因此， Tarlor 級數法在實際問題中很少應用，它的作用在于啟發我們去探索更好的方法。3-1 Runge-Kutta 方法的基本思想根據微分中值定理，存在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) 稱為區間xn,xn 1 上的平均斜率。由此可見，只要對這個平均斜率提供一種算法，由(3-2)式便可以得到一種計算公式。若取0，即取(xn, yn)點的斜率作為整個區間xn,xn 1上的平均斜率，這時K* f (xn, yn)，則得Euler 公式yn 1yn hf (xn , yn) n 0,1,2,若取1，即取(xn 1, yn 1)點的斜率作為整個區間xn,xn 1上的平均斜率，這時K* f (xn 1, yn 1)，則得后退的Euler 公式yn1yn 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) 兩點斜率的算術平均作為整個區間xn,xn 1上的平均斜率，則得下列公式，yn 1yn h (K1 K2)2K1f(xn,yn)( 3-3)K2f(xn1,yn hK1)這顯然是改進的Euler 公式(即梯形公式)。我們知道，改進的Euler 公式是具有2 階精度，而Euler 公式和后退的Euler 公式只有1 階精度。可見用兩個點斜率的算術平均作為平均斜率比只取一個點的斜率作為平均斜率精度要高些。由此得，如果設法在xn,xn 1內多預測幾個點的斜率值，然后將它們加權平均作為平均斜率K *，則有可能構造出具有更高精

15、度的計算公式。這就是Runge-Kutta 方法的基本思想。3-2 二階 Runge-Kutta 方法公式 ( 3-3) 是一個特殊的二階Runge-Kutta 方法。 它用了xn，xn 1這兩個點上的斜率值。我們現在在xn,xn 1內取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)其中含有三個待定參數1, 2, p。確定這三個待定參數的原則是使得公式 ( 3-5) 具有二階精度。為此將 K2在 (xn, yn )點作二元函數 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 級數法公式為h2yn 1 yn hynyn1, 2, p應滿足3-6)12112p 2滿足（3-6）形如（3-5）的所有公式均稱為二階Runge-Kutta 公1式。特別

17、當121 ， p=1 時就是改進的Euler 公式。3-3 三階 Runge-Kutta 方法為提高精度，我們現在在xn,xn 1內取xn，xnp， xn q0p q 1三個點，用這三個點上的斜率值K1 ， K2， K 3的線性組合作為平均斜率 K*。顯然，K1 f(xn,yn)，K2 f(xn p,ynphK1)，而K 3 的預測值我們用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）中的待定參數1, 2, 3, p,q,r,s。我們同樣用將K1， K2， K3作 Taylor 展開，然后代入（3-7）中的yn 1。通過與 三階 Taylor 級數法比較系數，則得待定參數滿足的條件為r s 1, 1231( 3-8)12p 3q 21以及2 p23q2( 3-9)23133pqs36參數滿足 ( 3-8)( 3-9) 的形如 ( 3-7) 的公式統稱為三階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 公式。下面就是一個經典的四階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 ) 。應注意的是，Runge-Kutta 公式的推導基于Taylor展開方法，所以它要求解函數有很好的光滑性。即y(x) 要具有所要求的導數 。 若 不 然 ， 用 高 階 Runge-Kutta

20、 公 式 可 能 不 如 用 低 階 Runge-Kutta 公式效果好。§ 4 單步法的收斂性與穩定性4-1 收斂性定義 2 若一種數值方法對于任意固定的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)有以下結論，定理 假設單步法有p 階精度， 且增量函數(x, y,h)關于y 滿足 Lipschitz 條件。即對任意y, y，總存在常數L>0, 使(x, y,h) (x, y,h) L y y( 4-3)則單步法(4-1 )收斂。且整體截斷誤差為y(xn) yn o(hp)。定理的證明從略。作為應用我們考慮梯形公式(改進的Euler方法)的收斂性。梯形公式的增量函數為(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)以上兩式相減，并注意到不等式的性質，得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)又由于我們已假設函數f(x,y) 關于y 滿足 Lipschitz 條件，所以上式變為(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所以，梯形公式的增量函數關于y 滿足 Lipschitz 條件。從而收斂。4-2 穩定性定義 3 若一種數值方法在節點xn上的值yn有大小的擾動，而以后節點xm (m>n)上的值ym所產生的擾動不超過，則稱 該方法是穩定

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

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

25、模型(4-4)變為預測：yn 1ynhyn校正：yn 1yn2h ynyn1h即yn 1 yn yn( ynhyn)2ynh2 ynyn2hyn1 221 h h y2n1yn* 11 h 12 2h2yn*1兩式相減得穩定性方程，n 1 1 h 1 2h2 n 得穩定性 條件為其他公式不再討論。例 見教材5 線性多步法線性多步法的公式為1 h 1 2h2 122 h0這一節我們研究線性多步法。解之得穩定區間為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）式中的待定參數，就可以得到一系列的線性多步法公式。確定待定參數時常用的方法有數值積分法和Taylor 展開法，在此我們用Taylor 展開法來 確定待定參數，因為這種方法更具有一般性，且更容易理解。設 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 項，所以應有下式成立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、個四步方法。先考慮顯式公式yn 10yn1yn 12yn 2h( 0yn 1yn 12yn 23yn 3)要使上式有四階精度，根據(5-2)式知其系數應滿足0121215-4)120123142214263131 8 2 3 112 2 27 31121231 16 2 4 1 32 2 108 3 17 個待定參數，而只有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)要使上式有四階精度，根據（5-2）式知其系數應滿足10125-6)1221012114221214211211218 23 13 112 211 16 2 4 1 4 1 32 2 17 個待定參數，而只有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 f ( x, y)例 解初值問題用顯式二步公式y(x0)y0yn 10yn1yn 1 h( 0yn1yn 1)。( 5-8)試確定待定參數0, 1，0, 1使方法階數盡可能高；并求局部截斷誤差。解 由于以上公式不易記憶，這里我們直接用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 級數法公式h2h3hp(p)yn 1ynhynyn ynyn2!3!p!為使方法階數盡量高，需有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 級數法公式當p=3 時相應的項為h yn o(h )4! n( 11)h4yn(4) o(h5)4!8作為練習請你用以上方法解下例：y f ( x, y)例 解初值問題用顯式二步公式y(x0 ) y0yn 1(yn yn 1) h( 0yn1yn 1)。試確定待定參數, ，0, 1使方法階數盡可能高；并求局部截