7常微分方程數(shù)值解法

1、第7章常微分方程數(shù)值解法7.0基本概念1. 一階常微分方程的初值問題y(x) = f(x,y(x), x(a,b)丿(7.0-1)y(a)= yo注：若f在D = a乞x zb , |y|+：內(nèi)連續(xù)，且滿足 Lip條件：L _0,使|f (x -yi) -f (x, y0| _L|yi y?|(7.0-2)則(7.0-1)的連續(xù)可微解y(x)在a, b上唯一存在。2. 初值問題的數(shù)值解稱(7.0-1)的解y(x)在節(jié)點Xi處的近似值yi ： y(xi)a X1 X2 . Xn = b.為其數(shù)值解，方法稱為數(shù)值方法。注： 考慮等距節(jié)點:Xi = a + ih, h = (b -a)/n. 從初始

2、條件y(a) = y0出發(fā)，依次逐個計算y1, y2,，yn的值，稱為步進法。兩種：單步法、多步法。 二階常微分方程y(x) = f (x, y(x), y(x)可設(shè)為一階常微分方程組的初值問題： 引進新的未知函數(shù) z(x) = y(x)，貝UJy(x)= z(x)z(x) = f(x,y(x),z(x)其初始條件為：:y(a) = y。Z(a) = y0稱為一階微分方程組的初值問題，方法類似。 邊界問題，常用差分方法解。7.1初值問題數(shù)值解法的構(gòu)造及其精度7.1.1 構(gòu)造方法對于(7.0-1)可借助Taylor展開(導(dǎo)數(shù)法)、差商法、積分法實現(xiàn)離散化來構(gòu)造求積公式:1. 設(shè) y Ca, b將

3、 y(xi+1)= y(xi+h)在 Xi 處展開y(Xi 1) = y(Xi) hy (Xi)FXi)hf(Xi,y(Xi)y()三Xi, Xi+1二 y(xi+i) ： yi+hf (xi, Yi)其中 yi ： y(xi).稱比+1 = yi + hf (Xi, yi). i = 0, 1, 2, ., n -1(7.1-1)為Euler求解公式，(Euler法)2. 用差商來表示:一 y(xi)茫yg y(xj).h得差分方程：出比二f(xi,yjh=yi+1 = yi + hf (xi，yi).即為 Euler 公式。若記 y(xi d:、y (Xi 1)= f (Xi 1, y(X

4、i 1)h=.yi+1 = yi + hf (Xi+1， yi+1).(7.1-2)稱為向后Euler法。注：Euler法為顯式，向后 Euler法為隱式須解出yi+1. 可用迭代法 yi+1(k+1) = yi + hf (xi+1, yi+1(k) k = 0 , 1 , 2, 解得 yi+1其中 yi+1 = yi + hf (Xi , yi).3. 對(7.0-1)兩邊取積分得x +yd) = y(Xi) + I f(x,y(x)dx(7.1-3)取不同的數(shù)值積分可得不同的求解公式，為：Xi + 用矩形公式：f(X, y(x)dx ： hf(Xi，y() ： hf 任 d,y(Xi J

5、)x=y(xi+1)： y(xi) + hf (xi, y(Xi) = Euler 公式y(tǒng)(Xi+1)： y(xj + hf (Xi+1, y(xi+1)= 向后 Euler 公式Xi -1hh2 用梯形公式：丨 f(x, y(x)dx= f (Xi, y(xj) + f (xr, y(Xi*) - f 牢，y)x212h=y(Xi1) ： y(Xi)- f(Xi,y(Xi)f(Xi1, y(Xi 1)2h= Y 1 ： Y - f (Xi, yi)f(Xi 1, yi 1)(7.1-4)2稱(7.1-4)為梯形公式一一隱式公式。弘=Yi +hf(x , Yi)顯化：預(yù)估值：hyr = Yi

6、十2【f (Xi,Yi)+ f (Xi+, Yi+)J厶校正值：.4. 幾何意義Euler法折線法改進Euler法平均斜率折線法例1:例 2:P473, P4747.1.2截斷誤差與代數(shù)精度定義7.1-1 稱1= y(Xi) -yi為數(shù)值解yi的(整體)截斷誤差。若yk= y(Xk) ,k = 0,1,2,，i-1.由求解公式得數(shù)值解y(Xi)，則稱ei=y(Xi)-為比的局部截斷誤差。注：局部截斷誤差是指單步計算產(chǎn)生的誤差，而(整體)截斷誤差則考慮到每步誤差對下一步的影響。定義7.1-2若求解公式的(整體)截斷誤差為O(h卩)則稱該方法是p階方法，或是p階精度。定理7.1-1設(shè)數(shù)值解公式：y

7、i+1= yi+ h(Xi,yi,h)中的函數(shù):(x,y, h)關(guān)于y滿足Lipschitz條件：| :(x, y,h) - (x, ,h)卜L | y - |，且其局部截斷誤差為hp+1階，則其(整體)截 斷誤差為hp階，即該數(shù)值解公式為 p階方式。注： 局部截斷誤差較易估計定理7.1-1表明：若ei = O(hp+1)則；i = O(hp).亠、h2.世2 Euler局部截斷誤差為ei 1y()=O(h)所以一階精度。2向后Euler法也是一階精度。 梯形公式為二階精度。例1:用Euler方法求解初值問題：y(x) = y(x) +(1 +x)y2(x), 1 ex v1.5.y(1) 1

8、1取步長h = 0.1，并與準(zhǔn)確解y( x)比較x解：因為 Xi = 1 + 0.1i，而 f(x, y) = y + (1 + x)y2，故2f(xi, yi) = yi + (2 + 0.1 i)yi于是Euler計算公式為yi+1 = yi + 0.1yi + (2 + 0.1i)yd , i = 0, 1, 2, 3, 4計算結(jié)果見P473表7.1-1注：Euler方法精度較低例2：用改進Euler方法求解初值問題：12y(x) = (y(x) y (x), 1 ex C1.5 xy(1) =0.5取步長h = 0.1,并與準(zhǔn)確解y(x) % 比較1 + x解：Xi = 1 + 0.1

9、 i,1 2f(Xi, yj (% - yi )二Xi(1 - yM1 -0.1i于是改進Euler法的計算公式為yiy0.1(1 -yjyi1 -0.1i+ 02 (1 -yjyi2 、1 -0.1i(1 - Vi 1)yi 11 -0.1(i +1) i = 0， 1, 2， 3, 4計算結(jié)果見 P474表7.1-2注：改進Euler方法精度比Euler方法精度高7.2 Runge Kutta 方法7.2.1構(gòu)造高階單步法的直接方法由Taylor公式：y(Xi J =y(Xj h)h 2h p.y(xi)hy(xi)-y(xi) . hp!y(p)(xi)p1y(p1)(p 1)!當(dāng)h充分

10、小時，略去Taylor公式余項，并以yi、yi+1分別代替y(xj、yg+1),得到差分方程:h2h p兒汀 hfgyjNfgyj.(p-1),(P (Xi.Yi)(7.2-1)其局部截斷誤差為:y(x J -i 1hp1(p 1)!(p 1)y即(7.2-1)為p階方式，上述方式稱為Taylor方式。注：利用Taylor公式構(gòu)造，不實用，高階導(dǎo)數(shù)f(i)不易計算。722 Runge Kutta 方法1.基本思想xi +因為y(Xi 1)= y(xj務(wù) f (x,y(x)dx=y(xi) + hf (，y()=y(xi) + hK ：其中K = f (，y()稱為y(x)在Xi，粕上的平均斜率

11、。若取 Ki = f (xi, yg)Euler 公式取 K2 = f (xi+i, y(Xi+i) 向后 Euler 公式一階精度1取一(Ki K2) 梯形公式二階精度2猜想：若能多預(yù)測幾個點的斜率，再取其加權(quán)平均作為K，可望得到較高精度的數(shù)值解,從而避免求f的高階導(dǎo)數(shù)。2. R K公式p%十 +匹qKjj 二 Ki = f(Xi, yj(7.2-4)jKj = f( +ajh, % +h送 bjsKs), j psm其中Kj為y = y(x)在 Xi + ajh (0 _aj_ 1)處的斜率預(yù)測值。aj，bjs，Cj為特定常數(shù)。3. 常數(shù)的確定確定的原則是使精度盡可能高。以二階為例：了卅=

12、yi +h(&K1 +沙2)K1 = f (Xi, yj(7.2-5)K2 = f(K +a2h, yi +b21hKJ(希望y(Xi+1)-yi+1 = O(hp)的階數(shù)p盡可能高)一方面：123y(Xi1)=y(Xi) hy(xj 二 h2y(xj O(h3)2!另一方面：將K2在(Xi, yi)處展開。(f (XoX yo：y)二 f(Xo, yo)(：xy ) f (x,y)excy2K2 = f (Xi, yi) + a2hfx(Xi, yi) + b21hK1 f y(Xi, yi) + O(h ).代入(7.2-5)得:23yi+1 = yi +hC1 f(Xi,yi)+hC2

13、f (Xi,yi)+ hC2a2 f x(Xi,yi) +b21K1 f y(Xi,yi)+ 0(h)23=yi + h(c1+ f (xi,yi) +C2a2hfx(xi,yi)+ (b21/a12)f (xi,yi)fy(xi,yi)+O( h )(希望)=yi h(C1 C2)y(x) C2a2h2y(xJ O(h3)希望：ei+1 = y(Xi+1)-yi+1 = O(h ).則應(yīng)：c1c2c?a 2b21=1-2=1a2特例：a2 = 1 - ci = C2 = 1/2 , b2i = 1，得 2 階 R-K 公式hyi十yi十2(心+心) 心=f(Xj, % )一心=f(Xi +h

14、, % +hKJ改進歐拉公式。C1 = 0 n C2 = 1, a2 = 1/2, b21 = 1/2，得：y = yi +hK2“ Q = f(Xi, yi)h hK2 = f(Xi +2, yi +2 Q)稱為中點公式。4.最常用的R-K公式一一標(biāo)準(zhǔn)4階R-K公式hy“ = yi + (K1 +2K2 +2K3 +K4)6(7.2-7)K1 =f (Xi,yi) K2 :hh二 f (Xi+c, yi+ K1)22hhK3 二=f (Xi+ ,yi+ -K2)22K4 :二 f (Xi+ h,yi+ ha)(7.2-8)輸入 a， b, n, yOh=(b-a)n，x0 = afor i

15、= 1, i=n, i+K1 = f(x0, y0)K2 = f(x0+h/2, y0+h*K1/2)K3 = f(x0+h/2, y0+h*K2/2)K4 = f(x0+h, y0+h*K3)x0 = x0+hy0 = y0 + h*(K1+2*K2+2*K3+K4)/6輸出x0，y0算法:function Runge_Kutta4(a,b,h,yO)n=(b-a)/h;xO=a;for i=1: nK1= f(xO,yO)K2=f(xO+h/2,yO+h*K1/2)K3=f(xO+h/2,yO+h*K2/2)K4=f(xO+h,yO+h*K3)xO=xO+hy仁 y0+h*(K1+2*K2

16、+2*K3+K4)/6;y0=y1en d;en d;fun ctio n f=f(x,y)f=x+y;en d;例1: (P478)用標(biāo)準(zhǔn)4階R-K公式求：(X)= x + y(x) 0 c x c 1y(o)=1的數(shù)值解。取h = 0.2，并與標(biāo)準(zhǔn)解y = 2ex x T比較。解：因為f(x, y) = x + y，從而由(7.2-8)得：yi +乎(心+2心+2心 *)K1+ Yi0.20.2*2=Xi+ + yi +K12丄0.2丄10.2K3=Xi+ + yi +K222K4=Xi+ 0.2 +y+ 0.2K3Xi+1K1K2K3K4yi+1y(Xi+1)0.211.21.221.4

17、441.2428001.2428060.41.4428001.6870801.7115081.9851021.5836361.5836490.61.9836362.2820002.3118362.6460032.0442132.0442380.82.6442133.0086343.0450763.4532282.6510422.6510821.03.4510423.8961463.9406574.4391733.436503注：步長h的選擇 (P479)使用數(shù)值解法求解初值問題 (7.0-1)，選擇步長h是一個重要問題。從每一步看，h小局部截斷誤差小，整體截斷誤差也就??；但從整個區(qū)間看，h小則

18、節(jié)點多，這不僅使計算工作量增大，而且也使舍入誤差的累積嚴(yán)重，導(dǎo)致數(shù)值不穩(wěn)定，所以步長h的大小要適當(dāng)。在實際計算中，如前所述，常常采用事后估計誤差及自動選擇步長的辦法來保證節(jié)點Xi處的數(shù)值解yi(i = 1 , 2,，N)滿足所要求的精度。設(shè)所用的一步法是 p階方法，其整體截斷誤差有漸近公式y(tǒng)(Xi+i) -yi+i = Mhp + O(hp+1)其中M是與h無關(guān)的常數(shù)。于是可用Richards。n外推法提高數(shù)值解的精度。事實上，從節(jié)點Xi出發(fā)，先以h為步長經(jīng)一步計算求出數(shù)值解yi(h1)，則整體截斷誤差如上式所示為yd)-y(h1 =Mhp +O(hp)(7.2-9)然后以h/2為步長，從節(jié)點

19、Xi出發(fā)經(jīng)過兩步計算求出數(shù)值解，則整體截斷誤差為(h) y(x 1)-yi21 =M+ 0(hp32 丿 用2p乘(7.2-10)式減去(7.2-9)式，得(h)(2p -1)丫以.1)-(221故有：(7.2-10)(h)p 1 - 1) =o(h )-1O(hp 1)(7.2-11)h(2)y(Xi+) yi 卅h(2)y 12p(h)yi 10(hp1)(7.2-12)1或:(h)y(Xi J 沖 12p2p -1(h)(y4 - y谿)+O(hp )(7.2-13)注: (1)(7.2-11)表明若取p (h)vyi 1 一(h)2卩尹，則整體截斷誤差為嚴(yán)，即精度提高階。(7.2-12

20、)表明數(shù)值解yt的整體截斷誤差近似為2p2-13)表明數(shù)值解丫畀的整體截斷誤差近似為y(h1)記二1_2p -1h(2)I yi +(h) |-yi 1 1，為精度要求，則y(Xi 1)h-yi +丹證玄芯且y(XiQ %(即步長h / 2是合適的；當(dāng)匚 ；時，說明步長h / 2仍然偏大，須將步長減半，繼續(xù)計算;當(dāng); ；/ 2p時，說明已有y(x屮)瘁2pb c s步長h / 2偏小，應(yīng)取步長h。例2(P481)用變步長的標(biāo)準(zhǔn)4階R-K方法求初值問題:y(x) =sh(0.5y(x) +x) +0.5y(x)o X v0.21.5y(o)=o的數(shù)值解，要求精度為;二1 10,27.3線性多步法

21、希望避免求多個點上f (x，y)的值，并且充分利用前面幾步的結(jié)果。一般形式：yi+i = ao yi + ai yi-i + +apyi-p + h(b-i yi+i +bo yi + + bpyi-p)(7.3-1)pp二 akyi上 h、bkyi丄 k=0k 一其中 yk= f (xk, yk), (k = i -p, i -p +1,，i, i+1)而 ak, bk為待定常數(shù)。 注：(1) 若b-i = 0時，(7.3-1)為顯式公式，否則為隱式公式。(2) 推導(dǎo)方法主要有：數(shù)值積分法和Taylor展開法。7.3.1數(shù)值積分法取節(jié)點：Xi, X/ , X Xi-3，作f的三次L-插值多項

22、式3P3(X)八k=03nj 式 Xi_jf (x* yi)，有 R3J4!3 n j j *記：f k = f(xk, yk) (k = i, i -1, i -2, i -3) xi-k = Xi -kh, x = Xi + th.代入上式有:3P3(X)二 p3(t)八k =0_ (t 1)(t - 2)(t3) f- 6 it(t - 1)(t2)-6仁所以:y(為 1)X 1T(X)x f(X，y(X)dXy(xjXi 15X P3(x)dx+O(h )y(xi)i50 p3(t)hdt O(h )t(t 2)( t 3) t(t 1)(t3)_fi 1 十_f-2 一 2y(xi)

23、h555f (Xi, yj -59f (Xi，yi)37 f (xg y(Xi “ - 9f(Xy(Xi)O(h ) 24(7.3-2)yi 1 二 y (55fi -59fi 137f 2 -9 3)24其中 f k = f(xk, yk) (k = i, i -1, i -2, i -3)稱(7.3-2)為4階Adams顯式公式。251h5f(5)()720若取節(jié)點：Xi+1，Xi，Xi-1，Xi-2，作f的3次L-插值多項式，可得 + hyi 1 = yi刃其局部截斷誤差ei彳=y(Xi -)i 4(9fi1 19fi -5fiA 仁)(7.3-3)4階Adams隱式公式:(7.3-4)

24、(7.3-5)其局部截斷誤差：e勺=y(xi -)-i - = - h5 f(5)()注： 隱式公式的顯化：(預(yù)測校正)Yi 1=yi一55fi -5937f/ -9仃 24h田彳化小“)19人一5仇f/24(7.3-6) 并非所有線性多步法公式(7.3-1)都可用數(shù)值積分法得到，但都可用Taylor展開法得到。7.3.2 Taylor 展開法設(shè)比十=y(xi -kh), yi-k = y(xi -kh)展開為:(-kh)j(j)yioOYu Tj壬j!j壬j!m 1(j) + (-k h)(m+) +(m 1)!J (j 1)Ti(j)yi 乂j =0J !(-kh)yi j# (j -1)

25、!代入(7.3-1)得：(kh)j(j)yi jm (j -1)!.(-如廣 y (m 1).yi+m!ak y +Zph bkk=AmzJ4m Lkh) j4 j!(j -1)!j(j)-yi(j)yipm h pakyiak( -k)Jk =0j 4 j! k =0p廣 bk(-k)jyi(j)k二hm 1pp- r ak(-k)m1、(m 1)! kk=為使(7.3-1)有m階精度，只須(7.3-7)的前m+1項與y(Xj+i)的展式:(7.3-7)y(N .J = y(Xi) hy (Xi)2!my (xjy(m)(Xi)O(hmd)m!對應(yīng)相等，即有方程組。pZ ak = 1k pp

26、、ak(-k)jj - bk(-k)j=1,(j =1,2, ,m)k k =此時有(7.3-8)hm 卅pp=y(x) -片=1-E ak(-k)m*-(m+1)E bk(-k)my(m*)(Xi) + O(hm42)(m+1)!k=0k=_1(7.3-9)特別取p = 3， m = 4有：a。+a*i *a2 *a3 =1a1 2a2 3a3 b + b。+ d + b2 + b3 = 13 +4a2 +9a3 +2b_1 2d -4b2 -66 =1-a1 -8a2 -27a3 +3b+36 +12b2 +27b3 = 1a1 16a2 81a3 4bj-4b -32b2-108b3 =

27、1令 ao = a1 = a2 = b-1 = 0，可彳得 a3 = 1 ， bo = 8/3， b1 = -4/3 ， b2 = 8/3 ， b3 = 0.代入(7.3-1)得Milne公式：4yi1=yiA 3h(2yi y2y)(7.3-10)3即4yi 1 =yi A ”h(2f (Xi, yj - f (Xy, yj 2f (x y)其中：e 1 二 y(Xi 1)- i 1 = 14 h5y(5)(Xi) o(h6)(7.3-11)45令 a1 = a3 = b2 = b3 = 0，得 Hamming 公式：13yi 1(9yiy h( % r 2yy)88(7.3-12)13(9

28、yi - y) h( f(Xi i, yi i)2 f(Xi, yi) - f (Xi，yi)881其中：e 1 二 y(Xi .J - i 1h5 y(5) (Xi) O(h4)3-13)40注： 與單步法相比，多步法不須反復(fù)計算f在(Xi, Xi+1)上某些點處的值，工作量大大減少。 多步法的前幾步須用同階單步法求之。例1用4階Adams顯式公式求初值問題：sh(0.5y+x)門廠 cc 廠y = +0.5y, 0 c x v 0.5j1.5y(o)=o的數(shù)值解，取h = 0.05。解：先用標(biāo)準(zhǔn) 4階Runge-Kutta公式求出此初值問題在x1 = 0.05， x2 = 0.1， x3

29、= 0.15處的數(shù)值解，然后用公式(7.3-2)求其余節(jié)點的數(shù)值解(p487) Milne 顯，Hamming 隱，顯化得 Milne-Hamming 公式：4%卅=y +(2fi -仁+2仁)3(7.3-14)13_yi 卅=(9% y) + h f (x,可卅)+ 2 人一f、一 8 87.4預(yù)估校正系統(tǒng)直接預(yù)估校正格式： 先用顯式公式算出預(yù)估值，再用同階隱式公式進行校正，沒有充分利用局部截斷誤差的信息。利用誤差補償?shù)霓k法，對預(yù)估值和校正值進行修正，可以使計算結(jié)果的精度更高一些。以4階Adams為例：251顯式：局部截斷誤差y(Xi1)-i1h5y(5)( )/(xi,Xi)(7.3-3)

30、隱式：局部截斷誤差y(Xj J彳=-9 h5y(5) ( ), - (XjX) (7.3-5)720設(shè)口+1, Ci+1分別表示Xi+1處數(shù)值解的預(yù)估值和校正值，則251 5 (5)丄 6y(Xi 1)- pi 1 = h y (Xi) O(h )(7.4-1)720y(Xi 1)y 1 - - h5y(5)(Xi) O(h6)(7.4-2)720兩式相減：cipi d = 270 h5y(5)(xi) O(h6)7205 (5)7206二 h y( J)(Ci 1 - 口 J O(h )(7.4-3)2702516代入(7.4-1)得 y(Xi 1) Pi 1 =270 (Ci 1 - Pi

31、 1)O(h )251 =y(Xi 1)-Pi 1(Ci 1 一 Pi 1)HO(h6)270即:251Pi 1(g 1 - Pi 1)作y(xi+1)的預(yù)估值，精度可提高一階。同理，將(7.4-3)代入(7.4-2)得196y(Xi 1)a 1 = (Ci 1 Pi 1) O(h )196=y(Xi 1) 一Ci 1 一 270 (Ci 1 一 Pi 1) =(h )(7.4-4)(745)19即：以Ci 1（c 1 - Pi 1）作為y（xi+1）的校正值，其精度可提高一階。270注：由于預(yù)估值ci+1還未算出，所以用前一步的q. Pi對Pi+1進行修正，由此得:hPi1 二yi(55fi

32、 -59仁 37仁 -9仁)預(yù)估242 51miPi 1 270 一 Pi)修正hCi 1 = yi (9 f (Xi 1, mi 1) 19 fi 5fi 4 fi )校正2419yi 1 二 g 1(g 1 Pi 1)修正270Pi 1 :=yi34h3 0i - f2仁)預(yù)估m(xù)i 1 :=pi 1112(G - Pi)修正8121Ci 1 :(9yi3h-yi)(f (Xi 1, mi 1) 2 fi -8fi j)預(yù)估yi 1 :=Ci1 _9 (、(Ci41 一 Pi 十)修正121注：在(7.4-6)和(7.4-7)中顯然無 C3，P3,可取 P3 = C3 = 0.(7.4-7)同理，對Milne公式（顯）和Hamming公式（隱）可得帶有誤差補償?shù)念A(yù)估校正系統(tǒng)公式:校正公式（7.4-6）求解例1:取步長h = 0.1，用帶誤差補償?shù)念A(yù)估y = y -紅，0 v x v0.5 yy（0） =1解：所需y1, V2, V3用標(biāo)準(zhǔn)4階R-K公式計算，然后按（7.4-6）計算g 5.P4927.5邊值問題的差分法基本思想：運用數(shù)值微分將導(dǎo)數(shù)用離散點上函數(shù)值表示， 從而將邊值問題的微分方程和邊界 條件轉(zhuǎn)化為只含有限個未知數(shù)的差分方程組， 并將此差分方程組的解作為該邊值問題的數(shù)值