科學計算與Matlabkj9

1、第九章常微分方程初邊值問題數值解同濟大學數學系計算數學教研室一階常微分方程（組）的初值冋題:(1)y = f (x, y),X e a, b,y(a) = yo,y(x)a, b上的m維函數向量，f(x,y)是定義m + 1維G = (x, y) | x e a, b, y e Rm上的m維已知函數向量.兩階常微分方程邊值問題:y + qy = f,x e a, by(a) = a, y(b) = 0其中q(x)和f(x)在區(qū)間a, b上連續(xù)，q(x) 0.一階常微分方程(組)初值問題：f (x, y)G ,y 滿(Lipschitz) ,x e a, b及y G Rm, y G Rm，總存在

2、L 0, 1|f(x,y) -f(x,y)| 0, 1|f(x,y) -f(x,y)| L|y - z|, HYPERLINK l bookmark2 o Current Document (1)存在唯一解，且解連續(xù)依賴于初始條件和右端.兩階常微分方程邊值問題： HYPERLINK l bookmark5 o Current Document (2)存在唯一解，且解連續(xù)依賴于邊界條件和右端.，其解y = y(x)都是區(qū)間a, b上關x.，其解y = y(x)都是區(qū)間a, b上關x.a = x0 x1 xN-1 xN = b 為.y (x)x” 貝yn.，其解y = y(x)都是區(qū)間a, b上關

3、x.a = x0 x1 xN-1 xN = b 為.y (x)x” 貝yn.,h” = x”+i x” = h = b-a，其h2.1歐拉公式及其改進-歐拉公式在節(jié)點X”處， HYPERLINK l bookmark2 o Current Document (1)式可寫成y (x”)= f(X”,y(x”).利用兩點數值微分公式y(tǒng)(x”) . y(叫 y(X”),可知y(x”+i) y(xn) + hf (xn, y(x”).y (x )是未知的，所以上式中y(xi)用其近似值yi代替.這就 得到歐拉公式y(tǒng)n+1 = y” + hf (x”, y”).(3)yo = y(xo),迭代可計算出各

4、節(jié)點處的近yi.算法2.1(歐拉方法)(1)N；h =強； y(1) = yo；k = 2,3,，N + 1做-x = a + (k 1)h;-y(k) = y(k 1) + hf(x, y(k 1);y.歐拉方法的M atlab程序如下：function x,y = odeeuler(f,yO,a,b,n) y(i) = y0;h = (b-a)/n;x = a:h:b;for i = 1:n,y(i+1) = y(i) + h*feval(f,x(i),y(i); end例2.1h = 0.1,用歐拉方法求解如下初值問題在區(qū)間0,1上的.y，= y - 2x，y (0) = 1.解:令yo

5、 = y(xo) = 1.由歐拉公式得到具體計算公式y(tǒng)n+1 = yn + h例2.1h = 0.10,1上的.(y = 1 + 2X)y = y - 2x， y (0) = 1.xn0.00.10.20.30.40.5yn1.00001.10001.19181.27741.35821.4351y (xn )1.00001.09541.18321.26491.34161.4142yn y (x”)|0.00000.00460.00860.01250.01660.0209Xn0.60.70.80.91.0yn1.50901.58031.64981.71781.7848y (xn )1.48321

6、.54921.61251.67331.7321yn y(Xn)|0.02570.03110.03730.04450.0527 HYPERLINK l bookmark2 o Current Document (1) y = y (x) xOy 平Po(xo, yo).；PoP0P1. Pn,拉折線法. HYPERLINK l bookmark2 o Current Document (1) y = y (x) xOy 平 HYPERLINK l bookmark2 o Current Document (1)y = y (x) xOy 平Po(xo, yo).；PoP0P1. Pn,拉折線法.

7、注意：誤差在不斷的積累和傳播！常微分方程yy = f (x, y)可改寫為dy(x) = f (x, y(x)dx.兩邊關于x在區(qū)間x”,x”+i上積分，得f X”+ly(x”+i) = y(x”)+f (x, y(x )dx xn要得到y(tǒng)(xn+1)的值，就必須計算右端的積分xn+1X”f (x, y (x )dx., . HYPERLINK l bookmark2 o Current Document 蔭方法.左矩形公式X”+lf(x, y(x)dx e hf(xn, y(x”).x”代入可得y(x”+l) e y(xn) + hf(xn, y(x”).y”y (x”)，y (x”+i)的

8、近y”+i，=yn+1 = yn + hf (xn, y”).右矩形公式xn+1xnf(x, y(x)dx e hf(xn+1,y(x”+i).可得到后退歐拉公式y(tǒng)n+1 = yn + hf(Xn+1, yn+1).xn+1xn梯形求積公式f(X, y(x)dx 2(f(Xn, y(Xn) + f (xn+1, y(x”+l).于是得到計算初值問題的梯形公式y(tǒng)n+1 = yn + 2(f (xn, yn) + f (xn+1, yn+1).(5)由于梯形求積公式比矩形求積公式的代數精度高，因此一般情況 HYPERLINK l bookmark9 o Current Document (5) H

9、YPERLINK l bookmark8 o Current Document (3).算法2.2(梯形方法)(1)N;h =強； y(1) = yo；k = 2,3,，N + 1做-x = a + (k 1)h;-y(k) = y(k 1) + 0.5h(f (x, y(k 1) + f(x + h, y(k);返回向量y.例2.2用梯形公式求解例2.1中的初值問題.: 2.1.:yn+1） = yn + 2仆-曽）+ 臨-愆）.:當h充分小時，該迭代法收斂于唯一不.例2.2用梯形公式求解例2.1中的初值問題.Xnyny (x”)yn y (xj0.01.00001.00000.00000.

10、11.09571.09540.00020.21.18361.18320.00040.31.26541.26490.00050.41.34231.34160.00070.51.41511.41420.00080.61.48431.48320.00100.71.55041.54920.00120.81.61391.61250.00150.91.67511.67330.00181.01.73411.73210.0021利用高階數值積分公式，還可構造其他離散格式.例如，阿達姆斯-巴士福德(Adams-Bashforth)公式y(tǒng)”+i =y” + 24(55f (xn，yn) 59f(Xn-1, y”_

11、i) +37f (x”_2, y”_2) 9f (x”_3, y”_3).阿達姆斯-莫爾頓(Adams-Moulton)公式y(tǒng)n+1:=yn + 24(9f(X”+1, yn+1) + 19f(xn, yn) 5f(Xn_l, y”_l) + f (Xn_2, y”_2)單步法與多步法，顯格式與隱格式-y”和y”+i,只要知道y”即可求y”+i，.單步法與多步法，顯格式與隱格式y(tǒng)”和y”+i,只要知道y”即可求y”+i，.-yn+1 時,y” ,yn-i, y”_2 等,:法.單步法與多步法，顯格式與隱格式y(tǒng)”和y”+i,只要知道y”即可求y”+i，.-yn+1 時,y” ,yn-i, y”_

12、2 等,:法.-yn+1,把右端yn+1,.單步法與多步法，顯格式與隱格式y(tǒng)”和y”+i,只要知道y”即可求y”+i，.-yn+1 時,y” ,yn-i, y”_2 等,:法.-yn+1,把右端yn+1,.-yn+1，求 y”+i 時,.f (x, y )關y,.局部截斷誤差可用于表征求解初值問題的數值方法的計算精度.一般地，常微分方程初值問題的數值解滿足形如yn+1 = yn + hg (y”+l, yn , yn-r),(6)，其中y”,，y”_r為y在r + 1個節(jié)點x”,，x”_r處的數值. HYPERLINK l bookmark10 o Current Document (6)，則

13、利用該式得到的數值yn+1y(Xn+1)”+1 = y(x”+l) y(xn) hg(y(xn+1), y(xn), , y(xn_r),被稱為數值方法的定義2.1若初值問題離散格式的局部截斷誤差為”+1 = O(hp+1),則稱該 離散格式公式所代表的數值方法具有P階精度或稱為p階方法.例2.3.:設y(x) a, b上充分光滑，M = max |y(x)|.x0 x由泰勒展開可得h2y(xn+l) = y(xn) + hy (xn) + _y (E), xn E xn+1 假設歐拉公式中右端的y” = y(xn),則kn+i|=|y (x”+i) y (x”) hf (xn, y (xn

14、)1=|y(Xn+1) y(X”) hy(Xn)|=h ly (E)I M h2 = O (h2),定義2.1若初值問題離散格式的局部截斷誤差為”+i = O(hp+1),則稱該 離散格式公式所代表的數值方法具有P階精度或稱為p階方法.例2.3.:設y(x) a, b上充分光滑，M = max |y(x)|.x0 x假設歐拉公式中右端的y” = y(xn),n+i Mh2 = 0(h2).這表明歐拉公式的局部截斷誤差是h2的同階無窮小，即歐拉公式.定義2.1若初值問題離散格式的局部截斷誤差為”+i = O(hp+1),則稱該 離散格式公式所代表的數值方法具有P階精度或稱為p階方法.例2.3歐拉

15、公式是一階方法.,-巴士福德公-2.1 2.2 , :梯形方法相對于歐拉方法雖然精度提高了，但是每一步都要用迭.,f (x)，丿.2.1 2.2 , :梯形方法相對于歐拉方法雖然精度提高了，但是每一步都要用迭 TOC o 1-5 h z .,f (x)，丿.,:h ，:. .yn+1 = yn + hf (xn, y”)，仃、yn+1 = yn + 2(f(xn, yn) + f (xn+1, yn+1),這就是一個預估校正公式. HYPERLINK l bookmark11 o Current Document (7)式通常稱為以下是改進歐拉方法的算法描述：算法2.3()(1)N;h =強；

16、 y(1) = yo；k = 2,3,，N + 1做-x = a + (k 1)h;-yp = y(k 1) + hf(x,y(k 1);-yc = y(k 1) + hf(x + h,yp);-y(k) = 0.5(yc + yp);y.改進歐拉方法的M atlab程序如下：function x,y = odeIEuler(f,yO,a,b,n) y(i) = y0;h = (b-a)/n;x = a:h:b;for i=1:n,yp = y(i) + h*feval(f,x(i),y(i); yc = y(i) + h*feval(f,x(i+1),yp); y(i+1) = 0.5 *

17、(yp+yc);end例2.4用改進的歐拉公式求解例2.1中的初值問題.:h = 0.1,2.1其具體形式為yn+1=yn + h (y” 鄉(xiāng))yn+1=yn + (yn -) + (y”+1 -) 例2.4用改進的歐拉公式求解例2.1中的初值問題.Xnyny (x”)yn y (xj0.01.00001.00000.00000.11.09591.09540.00050.21.18411.18320.00090.31.26621.26490.00130.41.34341.34160.00170.51.41641.41420.00220.61.48601.48320.00270.71.55251

18、.54920.00330.81.61651.61250.00400.91.67821.67330.00481.01.73791.73210.0058簡單的預估校正格式都包含兩個計算公式：, ；.當然也可以構造包含多個計算公式的預估校正格式.簡單的預估校正格式都包含兩個計算公式：, ；.當然也可以構造包含多個計算公式的預估校正格式., .-,.因, ,；.下面的預估校正格式y(tǒng)n+1 = yn-1 + 2hf(xn, yn),yn+1 = yn + 2 (f (xn, yn) + f (xn+1, yn+1)要好一些.四階阿達姆斯預估校正系統(tǒng)yn+1 = y” + 24(55f (xn, y”)

19、59f (x”-i, y”_i)+37f(Xn_2, y”-2)- 9f(Xn_3, y”-3), y”+l = yn + 24(9f (xn+1, y”+l) + 19f(xn, yn)5f (xn-1, y”_l) + f(xn_2, y”_2).歐拉公式可改寫成y”+i =y” + ki,k1 = hf (xn, y”).用它計算yn+1，需計算一次f(x, y)的值.y” = y (xn),yn+iy(x”+i)的泰勒展開式的前兩,O (h2).預估校正公式可改寫為yn+1 =yn + 2 k1 + 2 k2,k1 = hf (xn, yn), k2=hf (Xn + h, yn +

20、ki).用它計算yn+1，需計算兩次f(x, y)的值.yn = y(Xn),yn+1y(x”+1)的泰勒展開式的前三,O (h3).龍格-庫塔公式的基本思想：f (x, y)，彳, ,:較，,.一般的，顯式-弼的形式為/ryn+1=yn + 刀 3ikii =12(a2+ 021fn)h? + 0(啟)將y (xn+1)在點x = x”處的泰勒展開y (x”+1)=y (xn + h) = y (x”) + y(x”)h + y(x” )h2 + 0 (h3)=y(x”) + fnh + 1(d” + df fn)h2 + 0(h3).逐項進行比較上面兩式，為使y(x”+i) y”+i =

21、O(h3),則有1 + 2 = 1,宀2。2 = 2,2021 = 2.4a2,021,502., ,-令 a2 = 1,1 =2 = 2, 021 = 1,則有yn+i =yn + 2 k1 + 2 k2,k1 = hf (xn, yn), k2=hf (xn + h, yn + ki).這就是前面的預估校正公式.為使y(x”+i) y”+i = O(h3),則有1 + 2 = 1,宀2。2 = 2,2021 = 2.4a2,021,502., ,-1 = 0,2 = 1, a2 = 2, 021 = 2,則得到另一個二-yn+i =yn + k2,k1 = hf (xn, yn),k2=h

22、f (xn + 2 h, y2 + 2 ki).該格式稱為.-庫塔公式. :-y”+i=y” + 1(ki + 4k2 + k3),ki = hf(xn, yn),I k2=hf (xn + 2h, y” + 1 ki),I k3=hf (xn + h, yn ki + 2k2).三階三段Heun公式+ +1_323y”+i=y”+ 4(ki + 3k3), k1 = hf(X”, yn), k2 = hf (Xn + k3 = hf (Xn +標準四階四段龍格-庫塔公式y(tǒng)n+1 =yn + +1(k1 + 2k2 + 2k3 + k4),k1 = hf (xn, yn),k2=hf(Xn +

23、 2h, yn + 2ki),k3 = hf(Xn + 2h, yn + 2k2),、k4=hf (Xn + h, yn + k3).四階四段Gill公式y(tǒng)”+i=y” + 6(ki + (2 2)k2 + (2 + /2)k3 + k4), k1 = hf (xn, yn ),k2 = hf (Xn + 2 h, yn + 2 ki),k3 = hf (Xn + 2 h, yn + 並2 1 ki + 2 -護 k2), 、k4 = hf (Xn + h, yn g2 k2 + 2+2 k3)., ,計算函數值的次數之間的關系并非等量增加的，它們之間的關系:計算函數值次數r1234567r

24、2 8階數p1234456r 2, - . , -.算法3.1(-庫塔方法)(1)N;h =片； y(1) = yo；k = 2,3,，N + 1做-x = a + (k 1)h;-ki = hf(x,y(k 1)；-k2 = hf(x + 0.5h,y(k 1) + 0.5ki);-k3 = hf(x + 0.5h,y(k 1) + 0.5k2;-k4 = hf(x + h, y(k 1) + k3);-y(k) = y(k 1) + (ki + k2 + k3 + k4)/6;y.例3.1- 2.1 :f(x, y) = y 2x, y(0) = 1, h = 0.1.f,力代入標準- ,

25、彳yn+1 =yn + 6(k1 + 2 2 + 2k3 + 2k + k4),ki=0.1k2=0.1 yn + kk3=0.1 yn + kk4=0.1 yn + k34(xn + 0.5h)2yn + k14(xn + 0.5h)2yn + k22(Xn + h)yn + k3 丿.S-98999018O9OSZ*1ZE9903EZIOl9-99O9fOS002Z9-l991?3Z9l609-0080ssisi72i9,i9E99173I9I809-9HIE0KE616KI39W61691ZO9-9292 0696S8fl乙乙乙憶909-9910309SISt7iyi899I3H17IT

26、O9-9/99106Z0i79li7*TWSWIKTf09-9110901161792*122161792109-6820096913E8111zZ9l28Tl709-9ZIW03I9171796OI99171796OIIO9-0000000000000*100000000!OO|(ux)X- UA(nx)Xuxux驚創(chuàng)卵繆級串廠乙矽期孝報纟藩野喙竽蔚甸呦甸款務因竿皆務囪嗥卑物犁竿皆務囪嗥卑定義4.1若某數值方法對任意固定節(jié)點Xn,當h T 0時，有yn T y (x” ),則稱.定義4.1若某數值方法對任意固定節(jié)點Xn,當h T 0時，有yn T y (x” ),則稱., :yn+1 = y

27、n + hg(xn, y”, h),其中g(Xn, yn, h)稱為定義4.1若某數值方法對任意固定節(jié)點Xn,當h T 0時，有yn T y (x” ),則稱., :yn+1 = yn + hg(xn, y”, h),g(Xn, yn, h).f且僅是關于Xn,yn和h的連續(xù)函數.例如歐拉公式 的增量函數為g(x, y, h) = f (x, y).改進的歐拉公式的增量函數為g(x, y, h) = 1(f(x, y) + f(x + h, y + hf(x, y).P階方法收斂性的方法，并指明了局.定理4.1P階單步法yn+1 = y” + hg(xn, y”, h)g(x, y, h)在a

28、 x b, x y +x, 0 h ho上連續(xù)，y滿足,g(x,yi, h) g(x,y, h)| 0為正常數，則由該單步法得到的y”+i滿足如下誤差估:ChP|y(X”+1) yn+l| 廠(e(b-a)L 1).(12)證:令eo = y(xo) yo = 0,則有e”+i=y (x”+i) y”+i=e”+i + y(xn) + hg(xn, y(x”), h) y” hg(x”, y”, h)=e”+i + e” + h(g(xn, y(x”), h) g(x”, y”, h),e”+i = y(x”+i) (y(x”)+ hg(x”, y(x”), h)表示局部截斷誤.利用三角不等式

29、，可得p階單步法的整體截斷誤差滿足如下估計:注意到(b - a)L/1、(b-a)L(1 + hL)n+1 (1 + hL)hL_ 冬 (1 + hL)hL 1時，;A勺=O(hp+l)時，其整體截斷誤差 e” = O (hp).由于整體截斷誤差比局部截斷誤差低一階，因此在構造高精度的, .:當p 1時，;A勺=O(hp+l)時，其整體截斷誤差 e” = O (hp).由于整體截斷誤差比局部截斷誤差低一階，因此在構造高精度的, .,g(x, y, h) = f(x,y), jf (x, y),.-,f (x, y)滿足李普希茨條件，-g(x, y, h) = 1(ki + 2k2 + 2k?

30、+ k4).設李普希茨常數為L.因為ki = f(x, y),所以|ki(x,yi, h) ki(x,y2, h)|L|yi y2|,|k2(x, yi, h) k2(x, y2, h)|L|yi y2 + 號 ki(x, yi, h) #ki(x, y2, h)| W(1 + 2hL)|yi y21,|k3(x,yi, h) k3(x,y2, h)|L(1 + 1 hL + 1(hL)2)|yi y2|,|k4(x,yi, h) k4(x,y2, h)|L(1 + hL + *(hL)2 + #(hL)3)|yi y2.-,f (x, y)滿足李普希茨條件，-g(x, y, h) = 1(k

31、i + 2k2 + 2k? + k4).故g(x,yi, h) g(x,y2, h)| n)產Sm(m n),|Sn|,，I Sm1 W |Sn|,m = n + 1, n + 2, ，則稱該數值方法是絕對穩(wěn)定的.數值穩(wěn)定性的分析是相當復雜的，且總跟微分方程的右 f (x, y)，丿,-.特別地，一個約定俗成的方法是對如下的“模型方程”y=入y, (Re入 0)來討論數值方法的絕對穩(wěn)定性.歐拉公式的穩(wěn)定性：將歐拉公式用于模型方程可得y”+i = (1 + 入 h)y” = (1 + z)y”,z =入h.y”時有擾動6n，由必的傳播造成y”+i產生擾動值d”+i, J6”+i = (1 + z

32、 )6”,故要歐拉公式穩(wěn)定，只要I1 + z | W 1.,z = -1為圓心，1為.4.2收斂性與穩(wěn)定性-單步法的穩(wěn)定性后退歐拉公式的穩(wěn)定性：將后退歐拉公式用于模型方程可得y”+i = y” + 入 hy”+i,即1y”+i =y”.1 z于是有6”+i =6”.1 Z為使后退歐拉公式穩(wěn)定，只要111 ZI|z 1| 1.，丿Z = 1, 1.，.標準四階四段龍格庫塔公式的絕對穩(wěn)定區(qū)域為234Z 2 Z 3 z 4 + Z + + + + + 2 + 6 +24若入 0為實數，則由上式可近似得到2.78 z 0.，當入，-0.55 z 0,-3 z 0.一般而言,前面介紹的方法大部分都可平行

33、地推廣到方程組的情況. u 一 01u一 0v=-x 2xv+x +1例如，考慮如下一階常微分方程組x 6 0,10,初值條件為u(0)v (0)0,10進行N等分，記h = 10 .將歐拉公式應用于該問題, 則有Un+1Vn+1UnVn0-x21-xnUnVn0 xn + 1U0V0下面討論應用歐拉方法求解上述常微分方程組的絕對穩(wěn)定性. 簡單計算可知伏,n+1 一悅n+1）,必+1）T和6們,必）卩分別表示第n和n + 1步數值解 的擾動.由于矩陣-I -x 的特征值為入2 = e2nix，故可以證明 當|h入計 1（或h 0.1）時，（專）|，即該方法是絕對穩(wěn)定的.對于高階常微分方程初值問

34、題Jy(n)=f (x,y (x),y z(x ),y (n-1)(x),x Ga,b,y(a) = yo, y/(a) = y0，y(n-1)(a) = y0n-1),若記U = y(x), y(x),，y(ni)T,yo = %, y0，y0n-1) T,則可將高階微分方程化為一階微分方程組，從而利用前面介紹的.例如，考慮初值問題y = xy x 2y + x + 1, x G 0,10, y (0) = a, 0(0) = b.令y = u和y = v,則該初值問題轉化成一階常微分方程組 u 一 01u+一 0vx 2xvx +1x G 0,10,初值條件為u(0)v (0)定義5.1

35、HYPERLINK l bookmark2 o Current Document 稱為是剛性方程組，若對x e a, b,威立: (1) Re (“j (t) 0, j = 1, , m;max1.1j mmin1j m其中“ (t)是雅可比矩陣(f.s稱為剛性比.mx m例如，考慮常微分方程組 u 924 _u124 51v+初值條件為u1-15 cosx 專 sin x,、9 +l i , x G (0, 0 9 cosx + 3 sin xu(0)v (0)924這就是一個剛性方程組，其雅可比矩陣為為24 ，剛性24 51比s = |著=13.-u = 2e-3x e-39x + * c

36、osx, v = e-3x + 2e-39x | cosx.-.I, h |“刈冬 2.78,兇.h 2398 沁 0.0713.h = 0.05和h = 0.1,用標準四階龍格-庫塔公式求解該剛性問題.h=0.05xu(x)v (x)UnVn0.11.793061-1.0320011.712219-0.87031520.21.423901-0.87468091.414070-0.85501480.31.131575-0.72499841.1305230.72289100.40.9094086-0.60821410.9092763-0.60794750.50.7387877-0.5156575

37、9.7387506-0.51558100.60.6057094-0.44041080.6056833-0.44035580.70.4998603-0.37740380.4998361-0.37735400.80.4136714-0.32295350.4136490-0.32290780.90.3416143-0.27440880.3415939-0.27436731.00.2796748-0.22988770.2796568-0.2298511-.I, h |“刈冬 2.78,兇.h 2398 沁 0.0713.h = 0.05和h = 0.1,用標準四階龍格-庫塔公式求解該剛性問題.h=0.

38、1xu(x)v (x)UnVn0.11.793061-1.032001-2.6451697.8445270.21.423901-0.8746809-18.4515838.876310.31.131575-0.7249984-87.47221176.48280.40.9094086-0.6082141-934.0722789.35400.50.7387877-0.5156575-1760.0163520.000.60.6057094-0.4404108-7848.55015697.840.70.4998603-0.3774038-34989.6369979.870.80.4136714-0.32

39、29535-155979.4311959.50.90.3416143-0.2744088-695332.013906641.00.2796748-0.2298877-30996716199352: ,-. , ,. ，:, .定義5.2如果一個數值方法的絕對穩(wěn)定域包含z = Xh的整個左半平面,Re(z) 0,A- .,A-, J取，: ,-. , ,. ，:, .定義5.2如果一個數值方法的絕對穩(wěn)定域包含z = Xh的整個左半平面, TOC o 1-5 h z Re(z) 0,A- .,A-, J, .:-A-穩(wěn)定的； A-2;, -.,A-穩(wěn)定的.一段二階隱式龍格-庫塔方法：y”+i=y” + hki,k1 = f (xn + h, yn + h kl).二段二階隱式龍格-庫塔方法：y”+i=y”+ h(ki + k2),k1 = f(xn, yn),k2 = f (xn + h, yn + 2(k1 + k2).二段四階隱式龍格-庫塔方法:yn+1 =yn + h(k1 + k2),k1 = f (xn +(1 + g3)h, yn + h(k1 + (1 +k2),k2 = f (