計(jì)算方法常微分方程數(shù)值方法

1、第六章常微分方程數(shù)值方法連續(xù)問(wèn)題的離散處理一一尋求微分方程的解y(t)在某些離散點(diǎn)上的值一一y(t)在ti處的近似值yi ；記號(hào)：y(tj 所求函數(shù)在ii處的(準(zhǔn)確)函數(shù)值，yi計(jì)算得到的ii處的函數(shù)(近似)值，ti 節(jié)點(diǎn)，hi =1 -一一步長(zhǎng)；b 一 a等步長(zhǎng)：h三hi;以下若不另作說(shuō)明，一般總記h為等步長(zhǎng)n§6.1初值問(wèn)題的數(shù)值方法考察微分方程初值問(wèn)題：:y"(t) = f(t,y)a 蘭t 蘭b(、 (6.1) (a)=yo由微分方程理論可知：若函數(shù)f(t,y(t)關(guān)于y滿(mǎn)足Lipschtz條件，即存在與y無(wú)關(guān)的常數(shù)L，使f (t,yj f(t,y2)| 蘭L %

2、y21 Wt"a,b, yi,yc,d初值問(wèn)題的解存在且唯一；6.1.1 Euler法及其變形1、Euler 法由Taylor展開(kāi)式：h2h2y(ti 1) = y(ti) h y (ti) y ( i)二 y(ti) - hf(ti,y(ti) y ( i)22作局部化假設(shè)：yi=y(tj，并略去O(h2)項(xiàng)，便有y(ti 1Y y hf(ti,yi)將此式右端作為y(tj 1)的近似彳，便得到公式y(tǒng)i 1 = % h f (ti, yi) Euler 公式兩者的差(即略去的0(h2)項(xiàng))稱(chēng)為局部截?cái)嗾`差，記作E(ti, h):E(ti,h)二牛y ( i)2可有:Euler公式也

3、可由其他方法導(dǎo)出，例如由第四章數(shù)值導(dǎo)數(shù)公式,y(tiHy(ti1Vy(ti)hi ti ,ti 1 ；4 / 24解出y(ti 1)，并由y (tj = f(ti ,y(tj)替換，便可得Euler公式又如，根據(jù)Newt on-leib niz公式: y(ti J -y(ti) = / y (t)dtLi將y(t)以“ 0次插值多項(xiàng)式”替換，即以y (t) = y(tj y ( i)(t -tj代入積分，得到數(shù)值積分(左矩形)公式，及其誤差: h2t y(t)dt 二hy (tjy ( i)ti2又得到與Taylor展開(kāi)式相同的表達(dá)式，從而又導(dǎo)出 Euler公式。幾何意義E u l e折線法若

4、一個(gè)公式的局部截?cái)嗾`差為O(hp1)，則稱(chēng)該公式的精度為p階，或該公 式為p階公式。Euler公式是1階公式若取數(shù)值導(dǎo)數(shù)公式：y(tii)-y(ti) h 、 y(ti i)- y ( i)h2與前相同的推導(dǎo)過(guò)程，可以得到注意，以上的截?cái)嗾`差是在局部化假設(shè)的前提下得到的，即認(rèn)定yy(ti)倘若在每一步都按局部化假設(shè)，我們有Euler公式的總體截?cái)嗾`差：n h2bb - a巴E(y(b),h) = y(b)訃 八，y ( i)二hy ()y 222、后退Euler法h2”尹y(ti .J = y(ti) hf (ti i,y(ti i) y ( i)2在局部化假設(shè)的前提下截去局部截?cái)嗾`差h2E(

5、h,h) 一丁 y ( J2便得到后退Euler法公式:YiYi hf(ti i,yi i)注意到此公式中的右端也有yi i，需要求解關(guān)于yi 1的方程才能得到y(tǒng)i .1。因 此將這類(lèi)公式稱(chēng)為隱式公式，而將可以通過(guò)直接計(jì)算得到的公式稱(chēng)為顯式公式。后退Euler公式是一階隱式 公式，Euler公式是1階顯式公式6.1.2多步法1、梯形公式：ti +在式y(tǒng)C-y(ty(t)dt右端的積分中，取梯形積分公式，有h|h刑崑y(ti i) -y(ti) = y (ti i) y (ti) - “ y ( i)2 12由此，并據(jù)微分方程，可得：h梯形公式y(tǒng)i .1 = yi 2f(ti,yi) f (ti

6、 i, y .Jh3局部截?cái)嗾`差：E(ti,h) = 一y ( J12這是一個(gè)2階隱式公式。2、Simpson公式：t丄在式y(tǒng)(tj- y(tj二)y (t)dt右端的積分中，取Simpson積分公式，有% 1h 一h5y(ti i) y(tid) =：y(ti 1) 4y(tj y(ti)一£ y(i)390由此，并據(jù)微分方程，可得:Simpson 公式hyi 1 =yi4 - -f (tij,Yi j) 4f (ti,yi f(ti i, yi 1)3局部截?cái)嗾`差：h 5E(ti,hH-_y(5)( i)；90與以前的公式不同，用Simpson公式計(jì)算yi 1，必須有前2步的函數(shù)

7、值：yi 和yi。因此這種方法稱(chēng)為2步方法，而為啟動(dòng)此算法所需的最初的 2個(gè)函數(shù)值： yo,yi稱(chēng)為表頭。更一般的，若計(jì)算yi必須有至少前2步函數(shù)值，則這種方法 稱(chēng)為多步法。具體地，若計(jì)算yi，i必須有前k步的函數(shù)值，則這種方法稱(chēng)為 k步 方法,而為啟動(dòng)該方法所需的最初的k個(gè)函數(shù)值：y°,yi,，稱(chēng)為該方法的表 頭。與此相對(duì)，以前的方法計(jì)算yi，只須前1步的函數(shù)值yi，便稱(chēng)為單步方法。因此，Simpson公式是2步、4階、隱式方法。3、Adams方法(線性多步法)1 A在式y(tǒng)(t) - y(t=J y(t)dt右端的積分中，若取具有k+i節(jié)點(diǎn)的插值多 ti項(xiàng)式近似替代y(t)作為被積

8、函數(shù)，導(dǎo)出初值問(wèn)題的求解方法稱(chēng)為 Adams方法。(i)顯式 Adams方法Adams-Bashforth公式取tj, j二i,i -i, ,i -k處的yj二f(tjy)構(gòu)造插值多項(xiàng)式取代y (t):y (t) = p(t) R(t)iiP(t)八Tj(t)y (tj)Tj(t)f(tj,yj)j =s-kj 九-kR(t)二1(k 1)!y(k i 2)( iy (t)由于(t)-0,tti,ti i，有y(t G - y(l)=廣y (t)dt =廣E i j (t) f )dt + 廣yZ(W (t)dt匕 j 斗ti (k+1)!=送f +lj(t)dtf(tj,yjy(“中i) J

9、t%(t)dtj4± i(k 1)!J由此（以后為方便計(jì)，記 £ = f化，yj ）,可得顯式的多步法 Adams-Bash forth 公式及其局部截?cái)嗾`差yi 1 二 yihA(b0fi bi 仁bk fi _k )E(ti,h)二 rhk ：2 (k 2)y(i)8 / 24這是k 1步、k 1階的顯式公式F 表是 k =0,1，4 的 A,bj（j =0,1,，k）,r 的數(shù)值（注意：bj 二 A）kAb1b2bsb4r01112123-151221223-1653832455-5937-9251/72047201901-27742616-127425195288以

10、下是k =2的公式推導(dǎo)過(guò)程：(t 7 4)(t tj Q)(t ti)(t ti/)(t 十)(t -ti)P(t)f ifif i 2(ti t(ti 一t®(ty tj(ti4ti)(ti/ti)(ti<tj作變量代換：ti sh，以s為變量，當(dāng)t的變化區(qū)間為ti, ti 1時(shí)，s的變化區(qū) 間為0,1，且 dt = hds，有ti 1t P(t)dtti(S 1)(s2)2fi -s(s 2) fis(s 1)2fields(23fi-16f5fiJh4 1=y(4) ( i) 0s(s 1)(s 2)ds 3!(2)隱式 Adams 方法Adams-Moulton 公式若

11、取tj, j =i 1,i, i 一1,，i k 1處的= f(tj,旳)構(gòu)造插值多項(xiàng)式取代y(t),與前一樣的方法，可得隱式的多步法Adams-Moulton公式及其局部截?cái)嗾`差hyi 1 = yi A (b0 f (ti 1 , yi 1) 6 f i bk f i_k 1 )AE(h) "hk2y(k 2)( i)這是k步、k 1階的隱式公式。F表是 k =0,1，4的 A,bj(j =0,1,，k),r” 的數(shù)值(注意：b> A)：kAb0b；b2bTb4嚴(yán)011-1/21211-11221258-1-1124324919-51-17204720251646-26410

12、6-19-3160b述評(píng)：從表可見(jiàn)，對(duì)相同的k，A相同，而b：vbj,r*<r，特別是：<1， j jAIbJ(j =0,1，k)，而有若干j , - 1，因而在存在計(jì)算誤差時(shí)，由前步導(dǎo)致的誤差A(yù)顯然隱式公式要比顯式公式小(顯式公式對(duì)前步的誤差會(huì)被放大， 而顯隱式公式 則不會(huì))，而且局部截?cái)嗾`差也是隱式公式要比顯式公式小，結(jié)論： 隱式公式的 穩(wěn)定性一般比顯式公式好。6.1.3待定系數(shù)法利用Taylor展開(kāi)式比較有關(guān)項(xiàng)的系數(shù)，可以直接導(dǎo)出公式一一待定系數(shù)法。 例：求以下數(shù)值公式的系數(shù)r012,使公式具有盡可能高的精度：% 1 " % h(,fi十 J解：由于yi 1 y(t

13、i 1)，因此由Taylor展開(kāi)式，y(ti1)= y(ti h)二 yi hyi h2yi 卡“ h4yi O(h5)同時(shí)，對(duì)所求公式右端各項(xiàng)也作相應(yīng)的 Taylor展開(kāi)，并乘以相應(yīng)的系數(shù)：- y _ y - yih ：0£ 二二 h LoYih -1 fi j = h ：jy (tj -h)二 hh2 1- h3 y一-h4 -1 y(4) - O(h5)2 3!8心"侃心-2滬魂八注2嚴(yán)2宀仁4時(shí)(4)55)由于期望yi ：、y(tii)盡可能準(zhǔn)確，比較各對(duì)應(yīng)項(xiàng)的系數(shù)，可得方程組：E(ti,h)二 y(ti .J - yi -，因此注意到局部截?cái)嗾`差是E(h) h4y

14、i4!1 16 1 4 (4)8 5 , 4 (4)h y-h yi 3123! 120(h5) =3h4yj0(h5)8顯然，這就是k = 2的顯式 Adams-Bashforth公式。例：求以下數(shù)值公式的系數(shù)：“，、，空,使公式具有盡可能高的精度: yi 1 » y, h( 1?；?/)與上例完全一樣，比較系數(shù)，可得公式：4yi 1-h(2 f fi 4 2 仁)3局部截?cái)嗾`差：E(ti,h)二 14 h5y(5) O(h6)45這個(gè)公式稱(chēng)為Milne公式，稍后我們將用到它。綜上所述，從公式的構(gòu)造過(guò)程可見(jiàn)，微分方程初值問(wèn)題的算法公式基本上源 于：Taylor展開(kāi)式、數(shù)值積分或數(shù)值

15、微分，而插值方法一般也是通過(guò)數(shù)值積分或 數(shù)值微分完成的。Taylor展開(kāi)式待定系數(shù)法 直接來(lái)自于Taylor展開(kāi)式，通過(guò)比較系數(shù)，形成線性方程組，取得公式與局部截?cái)嗾`差；數(shù)值積分前面的多步法等，基本上都源于此，例如梯形 Simpsont i +方法等，而Adams方法用的是丫(如)=y(tj +y(t)dt，當(dāng)yt)用不同ti的點(diǎn)生成的不同的插值多項(xiàng)式代換時(shí)，便得到不同的公式，包括顯式和 隱式公式；數(shù)值微分由y(tj二也V y - y ( )導(dǎo)出Euler公式，而由h2y(tii)=川1八從) hy( i)可導(dǎo)出后退Euler公式，若取3點(diǎn)數(shù)h2值微分公式可導(dǎo)出更多公式：1h2小)訪®

16、;(t0) 4y(t1)-y(t2)l 亍()2h3- y(t 1H4y -3yd -2hf(tJ1yJr -y ()y(tlH2K- y(to)y(t2)-令 y ()h3=y(t 1HYj 2hf (t,yp y ()31hy(t2)埼 加-4曲1)3y(t2)l亍()3二 y(t 13l4y _yJ 2hf(t .1,y .1j29Ly ()3 9其中第二個(gè)公式就是“中點(diǎn)公式”；當(dāng)然它們也可以由 待定系數(shù)法取得：例 如第3個(gè)公式：用待定系數(shù)法求以下隱式公式的系數(shù)及局部截?cái)嗾`差：y 1 二 Ay ByChf (t 1, y .J 由Taylor展開(kāi)式：AAy(t)二By(ti J 二1-B

17、h3y.0(h4)6Chy Ch2y gch'yO(h4)y hy 1h2y - h3y - O(h4)2 6比較系數(shù)，得方程組：By_Bhy.1 Bh2 y.2Chf(t 1, y(t 1)二 Chy(t J 二比較：y(tJ =為使公式具有盡可能高的精度,A B =1_B C =1B2 _12a=43V/2:31 y(t 1) ： y 13與之相應(yīng)，有局部截?cái)嗾`差： 1 k>31 1 .31 2,3E(t ,h) h yh y61.6 32 3因此，公式是4y - y 42hf (ty J Lh3y0(h4)2 h3y O(h4)396.1.4問(wèn)題的性態(tài)與算法的穩(wěn)定性1、問(wèn)題

18、的性態(tài)一一問(wèn)題對(duì)原始數(shù)據(jù)的敏感性原始數(shù)據(jù)問(wèn)題的應(yīng)有的結(jié)果F（）原始數(shù)據(jù)擾動(dòng)問(wèn)題的結(jié)果F（）問(wèn)題的條件數(shù)相對(duì)誤差之比的上確界:_ m = Co nd (F)F(： )-F()/F)由微積分知識(shí):* M«St嚴(yán)）F0)Fg29 / 24初值問(wèn)題：固定步長(zhǎng)h，初值y0；計(jì)算結(jié)果：y（t）二F（y。）2 2F(yo) =y(t) =y(to h) =y(to) hy (to) O(h )hf(to,y°) O(h )因此:yo hf (to,yo)F (y。)： 1 - hfy(to, yo)Cond(F)=y(to)f氣 yo +hf(t°, yo)由于 y(to) ：

19、 y(toh)，當(dāng) h o -y(to)y(toh)(1 + hf；(to,y°)f；(to,yo? >o<o病態(tài)良態(tài)或：y(to) =yo擾動(dòng) (to) =yoy(t)二 yo hf (to,yo) O(h2)(t)二 (to) h (to) O(h2)二 (to) hf(to,y°、)O(h2)又 f(to, yo 、)= f(to,y°)fy(to,yo) 0(、2)得:二 (t) -y二(1 hfy(to，yo)因此：f(to, y°) £0 名£ 6良態(tài)fy(to,yo)0例:y(t)= 2y(t)方程的解例:y

20、(t)2y(t)方程的解例:y(t)= y(t)2y(t) 1-t T當(dāng)y <病態(tài)，> 6病態(tài)y（t）二ce2t 病態(tài) （圖示）y（t）二ce°t 良態(tài) （圖示）1fy =2y（t） -訂當(dāng)y ：丄 良態(tài)。（圖示）2t2、方法的穩(wěn)定性一一方法：由初始誤差導(dǎo)致的后期誤差的可控性。2 中點(diǎn)法：由數(shù)值導(dǎo)數(shù)公式八y(t hH )h y)2h6h3得 we “ r 2 hyty()可導(dǎo)出h 3中點(diǎn)法y =yi 二+2hf(ti,y，E(t, h)= 3 y) = O(h3)此為2步2階公式，但穩(wěn)定性差。例：八10y,y(0) =1真解：y(t)二 e J0t,y(1) : 0.45

21、3999 10 *取 2 種不同的表頭：yi=y(tj, y；2) hf(t0，y°)(即一步 Euler 法)，F(xiàn)表中列出對(duì)不同的步長(zhǎng)：1 n計(jì)算y(1)的近似值yN，yN0)按后退Euler法計(jì)算獲得，y N1)按中點(diǎn)法由初值及表頭y1獲得,yNT按中點(diǎn)法由初值及表頭y1 獲得，誤差放大則是比值(1) (2) y1- y1h后退Euler法yN0)(1) yN)(2) yN)誤差放大100.976562 審 103-0.266 *1034-0.238* 104574110 =_40.725657 101-0.179*102 0.594101190910；_40.477118*10

22、-0.179*10-0.566* 1001112110鼻0.456273 *10*0.619訊10鼻-0.546*10°10798說(shuō)明：NNE = E(y(b),h)八 E(ti.h)i ¥im1、總體誤差，由前可知后退 Euler法：y"G)=弓(b a)y“G) 丿 2(0,100, -t 0,1)可計(jì)算得E乞丄10九100=丄10，顯然實(shí)際 2 2因此，對(duì)于后退Euler法，若h=10,，(； (e'0t)E蘭丄10 *100 =丄；若h = 10*，可計(jì)算得2 2情況良好；而對(duì)于中點(diǎn)法，其總體誤差NNE 二 E(y(b),h)八 E(ti.h)八=

23、iimh3y'Vi) =2(b-a)y""G)h2<3考慮表頭y。,都用準(zhǔn)確值(注意在浮點(diǎn)數(shù)系中運(yùn)算，仍是有誤差的，這可認(rèn)為 是初始誤差)，若 h =10,L (e'0t)u(0,1000, X e 0,1)可計(jì)算得 |E <10/3， 而h=10,，可計(jì)算得|e卻0貽，而實(shí)際計(jì)算并非如此。2、對(duì)于同是中點(diǎn)法：由此表可見(jiàn)，兩個(gè)中點(diǎn)法的不同結(jié)果是由不同的表頭 導(dǎo)致的，實(shí)際的不同僅是y,，yf的不同，而它們之間的誤差(或不同)導(dǎo)致的y(1) 計(jì)算結(jié)果的誤差被放大5000倍甚至1萬(wàn)余倍。說(shuō)明中點(diǎn)法并不是好的方法，盡 管它是一個(gè)2階方法，但實(shí)際計(jì)算不如后

24、退Euler法這樣的1階方法。定義：k步方法，若存在常數(shù)K,h。，使-h乞h。y -V < K 0mjaxi yj- yj , i=k,k+1,其中y及y是按此方法分別由表頭y°,yi, 小及丑，,?kj計(jì)算得到，則 稱(chēng)此方法是穩(wěn)定的。此定義因未限制ho，故實(shí)際可要求ho > 0，因此本定義又一一“漸近穩(wěn)定性” 或“ 0穩(wěn)定性”。實(shí)際計(jì)算關(guān)心對(duì)確定步長(zhǎng)h， 一個(gè)算法，每步計(jì)算誤差是否被放大？無(wú)法討論一般方程研究對(duì)象：典型方程：y 二y(t) a 乞 t 乞 b, y(a)二 y° 其中：Re( ) ： 0。原因：對(duì)n維常微分方程組：Y (t)= F (t, Y

25、(t)，DyF (t,Y (t)為n階矩陣，它的特征 值反映了矩陣的性態(tài)，當(dāng)ReC ) < 0時(shí)，方程是良態(tài)，若Re( J _ 0時(shí)則是病態(tài) (工程中，通常稱(chēng)之為“穩(wěn)定性” ，Re( ) ： 0時(shí)，系統(tǒng)穩(wěn)定，否則，不穩(wěn)定)。定義：對(duì)于典型方程yW二y(t)，及確定的h二h，個(gè)k步方法，若由表 頭小及 乂肆1,皿匚計(jì)算得到y(tǒng)及y，存在常數(shù)c,0乞c ：1使 yi -yi 蘭C°m藥 yj _yj,75則稱(chēng)此方法對(duì)h是絕對(duì)穩(wěn)定的，全體h的集合稱(chēng)為該方法的絕對(duì)穩(wěn)定域。事實(shí)上，也可將此定義改寫(xiě)成Yi - Yi 蘭C1jX %與一y ,i=k,k+1,對(duì)于單步法則為：yi - I EC

26、丫口 -yj ,i =1,2,對(duì)典型方程,單步法總有yi T = p(' h)yi = p(h)yi例如：Euler 法：yi j = yih% =(1 h)yi， p(h) = 1 h1 - 1后退 Euler 法：丫“二'1, yi 1yi， p(h)=(1h)1 - hh(1 + 殆2)-2 + h梯形法：yiyi-(yi yi 1) , yi 1 =2 yi ， p(h - h ;2(1 -九 h J2 - h由于單步法y“-y“ = p（h）（yi yj因此，當(dāng)p（h） <1則方法絕對(duì)穩(wěn)定：Euler法：1+h|c1，在復(fù)平面上 是以-為中心，1為半徑的圓內(nèi)部；

27、 后退Euler法：一<1，即1 - h > 1，在復(fù)平面上 是以1為中心，1為1 -h半徑的圓的外部。為稍后介紹方便，對(duì)一般的k步方法，寫(xiě)成其中x' : jyij 1 -j fy 1=0j衛(wèi)j z0°yi 1 ry1 "川心k%上10仇1E u l e法:yi 1- yi- hfi = 0后退Euler法：yj1 - yi - h£ 1 = 0中點(diǎn)法：y1-yi - 2hfi = 0梯形法：yi 1 - yh /ri(fi 1fi02Simps on 法: yi 1h-yi-（fi 1 4ffi）二 0 ；kk即:例如i _k= 0定義一一算

28、法（“）的特征多項(xiàng)式：二（，h）h二（）k（ ） - j 2_0 k 機(jī) k 亠 一-：k/r.tkj=0k匚）=7 打2 = '0 .十；kj£kd分別稱(chēng)為算法（*）的第一、第二特征多項(xiàng)式。定理：對(duì)h，若算法（“）的特征多項(xiàng)式的全體零點(diǎn) 滿(mǎn)足i<1，貝S h使該算法絕對(duì)穩(wěn)定,i （i=1,2,，k）都全體這樣的h的集合，稱(chēng)為該算法的絕對(duì)穩(wěn)定域。例如：中點(diǎn)法的特征多項(xiàng)式根 2,2 = h ± v'1 +h2，由于h -0:11 h2 一1;可見(jiàn)，中點(diǎn)法總是不穩(wěn)定的。h ： 0:615預(yù)估-校正方法原因：顯式方法一一簡(jiǎn)單；但一般穩(wěn)定性不好；隱式方法求解復(fù)

29、雜；但一般穩(wěn)定性好;將兩者結(jié)合一一問(wèn)題：如何組合1、誤差估計(jì)一一兩個(gè)方法獲得同一y(t的不同近似yi，以此估計(jì)誤差(1) 同階、同步長(zhǎng)的兩個(gè)不同方法方法 a： y(ti J _ 射 i ： hp 1y(p ")( 1)方法b: y(tiQ-片二歸卩裁卩刊(J)由于步長(zhǎng)h 一般較小，因此可認(rèn)為 y(p1)( ih- y(p1)( 2)，(實(shí)際上就是估計(jì)此 y(p1)( i)的值，因?yàn)橐坏┇@得此估計(jì)值，便可得到E(t, h) = y(ti .J - y： i的估計(jì)值)將兩式相減，得% i - yi i(：- ' )hp dy(p1)()即因此hp1y(p1)(')yi 1

30、 - yi 1(:- )y(ti 1)- yi 1(2) 不同階、同步長(zhǎng)的兩個(gè)不同方法，設(shè)方法 a： y(ti J -% 1 =Jhp1y(P "(j方法 b： y(ti1)-yi1 二 My(q1)( 2)由于hq1是hp1的高階小量，所以?xún)烧呦嗉?減)時(shí)，后者可略去；將以上兩式 相減，可得P 1 (P 1) y yi 1 - yi i - h y ( 1)因此：y(ti J - yi 1 ： % 1 - yi 1(3)同階、不同步長(zhǎng)的同一方法，設(shè)步長(zhǎng) h : h步長(zhǎng)h :y(ti 1) - yi 1= hp 1y(p 1)( 1)步長(zhǎng) h : y(ti 1)- / 1 = ： h

31、p 1 y(p 歸(2)仿(1)，將兩式相減，得yi1"1 j1-(hh)p1hp1y(p1)(1)所以y(ti 1) - yi 1 : (yi 12、預(yù)估-校正法(Predictor-預(yù)估,Evaluation-計(jì)算,Correct-校正)(1) Heun 方法改進(jìn) Euler法預(yù)估：Euler 法一Pi .1 二 hf (ti ,yj校正：梯形法 f(ti，yi)f(tii，pii)保持二階精度：y(ti J y =0(h3)證明：記討為采用原始的梯形法獲得的結(jié)果，即：y(tii)-y”=O(h3)，由本節(jié)初所作的假設(shè)：函數(shù)f (t,y)關(guān)于變量y滿(mǎn)足Lipschitz條件，y(

32、ti Jy(ti +) iyi +f(ti,yi)十f(t+yJ+£f(ti和PiG f(ti和門(mén)打2 2h*h沖<y(ti 卅)-% +qf (yj + f(t“,y ),+-f(ti41, PR f(ty,y )9h3)ELPv -y<O(h3- Pr y(t) +|y(tid y用、0(h3)3、h(2) Adams-Bashforth-Moulton 方法取k =3的Adams-Bashforth顯式公式作預(yù)估，Adams-Moulton隱式公式校正: hPi 1 二 yi -9fi37 fi/ - 59 fi55 fi24hyi i 二 fi - 5fid 19

33、 fi 9f (ti 1, Pi i)24對(duì)應(yīng)局部截?cái)嗾`差公式：y(ti i) - Pi i =禦 h5y(5)( i),y(ti i) - yi i h5y2)，7207 2 0根據(jù)前已介紹的誤差估計(jì)方法，有19y(ti 1) 一 yi 1(yi 1 - Pi 1)。27 0由于在計(jì)算過(guò)程中，右端的值均已得到，故在每一步，均可對(duì)誤差進(jìn)行監(jiān) 控。在計(jì)算中，對(duì)于給定的相對(duì)誤差界 Re，19 |yi 卅 一 Pi41a Re270y“| + Sm19yi+ 一 Pi卅Re270|yi+ Sm100若t h e n2h= h此處 Re(Relative error)相對(duì)誤差界,Ae(Absolute

34、 error)絕對(duì)誤差界,Sm(二Small)是一個(gè)與A%e相當(dāng)?shù)男≈?，目的是避免yi 0時(shí)出現(xiàn)的估計(jì)失真(A%©實(shí)際上是函數(shù)y(t)的數(shù)量級(jí)，只須注意y-y 一 y若相對(duì)誤差過(guò)大，須步長(zhǎng)減半(h/2二h)，新的進(jìn)程中，將由yi,yi_,等原 有的值計(jì)算y(ti h2)的近似y 1，而根據(jù)本方法，它是fi,f 1,fi,f 3的組合,2iti i 2 2 2對(duì)此中的原未取得的f i與f 3的值，可采用插值法得到，但應(yīng)注意，本方法是4階方法，因此，用插值法取得這些數(shù)值時(shí)應(yīng)保證誤差也應(yīng)至少是0(h5)的，而 5 點(diǎn)插值多項(xiàng)式的余項(xiàng)： R(t)二 fto,ti ,t2,t3,t4,t(t

35、-to)(t -ti)(t -t2)(t -t3)(t - t4), 且t -tj =0(h)，因此用此5點(diǎn)插值方法可以實(shí)現(xiàn)誤差是 0(h5)的。為此，取最 近相鄰的5點(diǎn)構(gòu)造插值多項(xiàng)式，可得所求點(diǎn)處f (t,y)的近似值：1f 3(3£_4- 2 0_39fi0 2f6_S fi 5 )匕 1281f 1(-5f 28fiJ3 -70fj2 140£35fJ ；i-2128(3) Milne-Simpson 方法以Milne公式為預(yù)估，Simpson公式做校正。4Pi 1 =y-h2fi -i 2 fi fij 4 fi( t1, P1)他們的局部截?cái)嗾`差公式分別為:y(t

36、iPi h5y(5)( 1),45y(ti J yi 1 = -1h5y(5)( 2)；90與前相同，兩者相減有:可導(dǎo)出預(yù)估公式-Milne公式的誤差估計(jì)28與數(shù)值積分的Romberg公式的導(dǎo)出思想一樣，現(xiàn)在設(shè)法將這部分的誤差“歸還” 給預(yù)估的Pi!。由于在得到Pi 1后，希望對(duì)Pi 1進(jìn)行修正，使修正后的值更接近 y(ti J再進(jìn)入校正公式參與計(jì)算，這樣得到的校正值有望更符合 Simpson公式的 實(shí)際效果。但由于尚未計(jì)算yi 1，只能用yi - Pi近似替換yi 1 - Pi.(實(shí)際上， 是 y(5) ( i) : y(5)( i J)，即: y Pi % 1 - Pi 1，因此有4Pi+

37、 = yi_3 + §h2fi_fi_t2fi丄28m = 口十+二(yi Pi)29hyi=yi斗十才戸口十彳右+ f (ti41,mr1)注：這是4步方法，必須有表頭y°,yi,y2,y3才能開(kāi)始計(jì)算。預(yù)估步與校正 步當(dāng)i =3便可執(zhí)行，但在第一次的修正步計(jì)算 m4中，由于此前無(wú)p3，因此只能 令m4二P4，也即第一次因?yàn)闊o(wú)法修正所以不進(jìn)行修正。(4)修正 Milne-Hamming 方法解初值問(wèn)題的Hamming公式：13y i (9y - yQh 2f, f (t, y, 1)88E(h,h) =y(ti i)-yi ih5y(5)()40以Milne公式為預(yù)估，H

38、amming公式進(jìn)行校正。與前相同，兩者的局部截?cái)嗾` 差相減有：1215(5)巴yi 1 Pi 1h y ( 1)360可導(dǎo)出預(yù)估公式-Milne公式的誤差估計(jì)1 1 2y(ti 1 ) - Pi 1 ：1 - Pi 1)1 2 1校正公式-Hamming公式的誤差估計(jì)：9y(ti 1)-yi 1(yi 1 - Pi 1)121與前一樣，現(xiàn)在設(shè)法將誤差“歸還”給預(yù)估的Pi 1，以及y1 ;由于通常將最終結(jié)果記作yi 1，因此將Hamming校正公式計(jì)算結(jié)果記作c, d ；與前Milne-Simpson 方法一樣，在修正Pi 1時(shí)，Ci 1尚未產(chǎn)生，因此只能使用Ci，從而形成公式：4Pm =+h

39、2fi - fi_ +2fi+ 112/ 、mi+ = P" + ,一(G 一 Pi)12113&十=匚(9% yi( +：h fi+ 2人 + f (t+m)889yi 1 =Ci 1 777(G 1 - Pi 1)L121注：仍與 Milne-Simpson方法一樣，這是一個(gè) 4步方法，必須有表頭 y°,y1,y2,y3才能開(kāi)始計(jì)算。預(yù)估步與校正步當(dāng)i =3便可執(zhí)行，但在第一次的修 正步計(jì)算m4中，由于此前無(wú)P3，因此只能令m4二P4。注：如果從求解隱式方程的迭代法來(lái)看，這校正步可以認(rèn)為是迭代了一步， 而迭代的初值則是預(yù)估值，或是預(yù)估的一個(gè)修正值。因此從迭代收斂

40、的角度看， 應(yīng)對(duì)初值提出一定的條件，通過(guò)研究，可以得到如下結(jié)果，各方法對(duì)步長(zhǎng)的限制:A-B-M4 方法：h ：0.75fy(t, y)Milne-Simpson 方法：_0.45h :fy (t,y)修正 Milne-Hamming 方法：-0.69h ：fy(t,y)6.1.6 Runge-Kutta 方法將Heun方法（改進(jìn)Euler法）改寫(xiě)成如下形式： 1y =yi +?（心十心）“心=hf （ti,yjK2 =h f（h +h,yi +KJ此為一單步、2階、顯式公式（比較梯形：2階隱式）。若注意此處的KlK2，它 們是兩個(gè)不同的增量，Heun方法則是從yi出發(fā)加上若干個(gè)增量的加權(quán)平均值

41、后 得到y(tǒng)i 1，將此推廣到一般形式，y出發(fā)加上m個(gè)增量的加權(quán)平價(jià)取得y： “，便 有：m -級(jí) Runge-Kutta 公式y(tǒng)i 十=yi *lKi +-2«2 十+-mKmKi =hf（ti,yj“K2 =hf（ti +Sh, yi +021心） <A A A A J A AKm =hf （ti *mh,yimiKi 行2心mmKm）為導(dǎo)出該公式，先引入 二元函數(shù)Taylor展開(kāi)式：Affc1cC nf(t h,y k)= f(t,y) (h k )f(t,y)k )nf(t,y)t: yn!t.y133(h k jn 1 f(t th,y =k), Q< 齊，出乞

42、1 (n 1)!；t:y以二級(jí)Runge-Kutta公式為例：J-Yi 1 *1心2心Kh f(ti,yi)K2 二 h f *： h, yiK1)先考慮準(zhǔn)確值y(ti J的展開(kāi)：123y(ti 1) =y(ti h) “ihyi h y O(h )2其中y = fY ' d f(t, Y(t)二 ftfyY = ft f fydt因此1 2 3y(ti .J =yi hfi 匕人(f/ f fy) O(h )另一方面，考慮近似值yi .1的展開(kāi)，由二元函數(shù) Taylor展開(kāi)式，二級(jí)Run ge-Kutta 公式中，K2 Mt ： hft0fy O(h2)?故應(yīng)有：=yihfi phf

43、j 2： h2ft 2：h2ffy 0(h3)將準(zhǔn)確值y(ti .J與近似值 1展開(kāi)式作比較，應(yīng)爲(wèi)+為=1二 1 k廬=_2若令2 = C，貝U'1 = 1 - C,變形Euler法:取 c，貝U 、二匕二、，=：=1 ,Heun 方法(改進(jìn) Euler 法)取c =1，貝卩、=0,，2=1, 二：二1?=yi +K2“ K1=hf(ti,yJh K1© = h f(ti +二+). 2 2(注意此公式與中點(diǎn)法的區(qū)別)以上兩公式均為二階二級(jí)Ru nge-Kutta公式。常用的四階四級(jí)Runge-Kutta公式(稱(chēng)為經(jīng)典、古典或標(biāo)準(zhǔn)Runge-Kutta公式Classical

44、Run ge-Kutta)1 y = yi +(Q +2心+2心十心)6心=hf(ti,yj丄 hx K1心=hf * ： %2 2丄h丄K2&=hf(ti + , yi +)2 2K4 二hf(ti h, yi K3)幾何意義：若將上述RK -4-4公式改寫(xiě)為：YiYi J6( f1 2f2 2f< f4)其中fj (j =1,2, 3,4)為分別為原式中Kj(j =1,2,3,4)中相應(yīng)的函數(shù)fj值。比較在區(qū)間tid上的數(shù)值積分Simpson公式ghyiyi - f (t, y)dt *(右4i6其中fi, fm, fi 1分別表示f (t,y)在區(qū)間ti,tii的左、中、右

45、端點(diǎn)的值；它們 分別有：fi:-fl,fm、ff3,fif4；因此可簡(jiǎn)略地得到該RK 一4一4公式的 誤差公式是E(h,h)誤差估計(jì)：在實(shí)際計(jì)算中，由于Runge-Kutta方法是一個(gè)高精度的單步方法。因此，通常作為求解常微分方程的一個(gè)獨(dú)立的方法使用。在計(jì)算時(shí)，常常要對(duì)計(jì)算值進(jìn)行誤差估計(jì)。以四階四級(jí)的標(biāo)準(zhǔn) Run ge-Kutta公式為例：若用同一公式的不同步長(zhǎng)方法進(jìn)行誤差估計(jì)，例如用步長(zhǎng)h和h2計(jì)算獲得y(ti 1)的不同近似進(jìn)行，為此必須計(jì)算11個(gè)f(t,y)值(其中：丫廠y1計(jì)算4一個(gè)比較簡(jiǎn)單的方法是Ru nge-Kutta-Fehlberg方法X),592K5K5055yi 1 二 y

46、i C25K121616yi 1 二 yi (K1135二h f(t ,yj1二 hf(tih, yi43=hf 魚(yú) 3 h, yi812二hf(tih,13K1K2K3K414082197K3K4 -25654104665628561K3k412825564301-K1)439K1K2)323219327200yiK1K21972197E(t,h) =0(h5)6),E(t,h) = O(h6)7296K3)21978454393680K1 -8心K3K4)2165134104835441859yi K1 2K2K3K47i 27 12 2565 3 4104 411282197 丄 1 丄

47、 _E(ti,h) ： y“ -丫K1K3K4K5心3604275752405055K5K6=h f (ti h, y1=h f (ti h,11Ks)402 一個(gè) f(t,y)值，yi y 1 yi 1 計(jì)算 8 個(gè) f (t,y)值，而 K1 =hf(ti ,yj 是重復(fù)的)。 2這比不估計(jì)誤差增加7個(gè)f (t, y)值的計(jì)算量以此作誤差估計(jì)，在計(jì)算過(guò)程中僅需計(jì)算 6個(gè)f (t, y)值的計(jì)算量。小結(jié)：常微分方程求解的理論問(wèn)題，主要是問(wèn)題的性態(tài)和方法的穩(wěn)定性。 算法包括 基本算法的構(gòu)造(或：形成)及組合。前者就是通過(guò)Taylor展開(kāi)式 Runge-Kutta 法、待定系數(shù)法形成如 Miln

48、e法等，或數(shù)值積分(例如：Adams方法，梯形法， Simps on法等)，或數(shù)值微分(或稱(chēng)數(shù)值導(dǎo)數(shù))形成；當(dāng)然同一方法實(shí)際上可以 通過(guò)多種方法導(dǎo)出。后者就是利用前面已得到的基本方法，進(jìn)行適當(dāng)組合，一方面利用顯式方法簡(jiǎn)單，作為預(yù)估，利用隱式方法比較穩(wěn)定，作為校正；進(jìn)一步， 同階的不同方法，或同一方法不同步長(zhǎng)的局部截?cái)嗾`差比較， 將誤差中的估計(jì)項(xiàng) y(p 1)()變換成可計(jì)算的 1 - 口 1 ,得到實(shí)際誤差估計(jì)值，或直接用來(lái)比 較誤差是否滿(mǎn)足要求，或用來(lái)對(duì)對(duì)公式作修正。6.1.7常微分方程組、高階常微分方程1、常微分方程組初值問(wèn)題：常微分方程組初值問(wèn)題的解法與常微分方程初值問(wèn)題的解法幾乎是一模

49、一 樣的，唯一區(qū)別就在于所有的公式都用向量代替純量。為此，記向量y2(t)費(fèi),F (t, Y(t) =5(1 y1 (t), y2(t),yn(t)" f2 (t, ydt), y2(t),yn(t)<yn(t)j<fn(t, y1(t), y2(t),yn (t)Y (t)二則可將常微分方程組初值問(wèn)題"(t) = f,t,y1(t), y2(t),yn(t) y2(t) = f2(t,y1(t), y2(t),yn(t)I- « « «yn(t) = fn(ty(t), y2(t),yn(t)y1(a)二 y10 , y2(a)二

50、 y20 , ,yn(a)二 yn0 ,改寫(xiě)成簡(jiǎn)單的向量形式：Y (t)二 F (t,Y (t), Y(a)二 Y°.對(duì)此向量形式的初值問(wèn)題的解法，只需將以前一維方程的數(shù)值方法改為向量形式 便可。例如，解常微分方程組的Euler方法，寫(xiě)成向量形式:Yi i 心 hF(ti，丫,標(biāo)準(zhǔn)Runge-Kutta公式的向量形式： 1y“ =yi + - (Ki + 2 K2 + 2K3 +K4)Ki = hF(t i ,Yi )h K1 © = h F(t-,Y-1)22hK2& = h F(ti +；,Y +嚴(yán))22I & = hF(ti + h，丫廣心)2、高階方

51、程：然后用對(duì)于高階常微分方程初值問(wèn)題，一般可以將它改寫(xiě)成常微分方程組, 前述的向量方法求解。例如：3階常微分方程初值問(wèn)題y 二 f (t, y(t), y (t), y (t), y(a) = y。，y (a) = y。，y (a) = y。，令V(t) =u(t)9(t)、u "(t) =v(t),Y (t) =u(t)vt) = f(t, y(t), u(t),v(t)lv(t)并記'u(t)“F (t, Y (t) =v(t),Y 0 =y0寸(t,y(t),u(t), v(t)丿<v0r>則原高階方程就可以改寫(xiě)成如下形式，從而可按標(biāo)準(zhǔn)方法求解，Y 二 F (t,Y (t), Y 二 Y。6.2邊值問(wèn)題數(shù)值方法簡(jiǎn)介6.2.1差分方法討論常微分方程：y (t) q(t)y(t)二 f(t) a t 汕，其中q(t) 一0 ;對(duì)于一般方程：y (t) p(t)y (t) q(t)y(t) = f (t) a 罕 5可以通過(guò)以下方法轉(zhuǎn)化為上述形式。對(duì)方程乘以e p(t)dt，得: