龍格庫塔方法RungeKutta

《龍格庫塔方法RungeKutta》由會員分享，可在線閱讀，更多相關(guān)《龍格庫塔方法RungeKutta（19頁珍藏版）》請?jiān)谘b配圖網(wǎng)上搜索。

1、龍格-庫塔方法（Runge-Kutta） 3.2 Runge-Kutta法 3.2.1 顯式 Runge-Kutta法的一般形式 上節(jié)已給出與初值問題(1.2.1)等價(jià)的積分形式 　　　　　　　　　 　(3.2.1) 只要對右端積分用不同的數(shù)值求積公式近似就可得到不同的求解初值問題(1.2.1)的數(shù)值方法，若用顯式單步法 　　　　　　 　　(3.2.2) 當(dāng)，即數(shù)值求積用左矩形公式，它就是Euler法(3.1.2)，方法只有一階精度，若取 　　(3.2.3) 就是改進(jìn)Euler法，這時(shí)數(shù)值求積公式是梯形公式的一種近似，計(jì)算時(shí)要用二個(gè)右端函數(shù)

2、f的值，但方法是二階精度的.若要得到更高階的公式，則求積分時(shí)必須用更多的f值，根據(jù)數(shù)值積分公式，可將(3.2.1)右端積分表示為 注意，右端f中還不能直接得到，需要像改進(jìn)Euler法(3.1.11)一樣，用前面已算得的f值表示為(3.2.3)，一般情況可將(3.2.2)的表示為 　　　　　　　　　　　　　 　　(3.2.4) 其中　　　　　 這里均為待定常數(shù)，公式(3.2.2)，(3.2.4)稱為r級的顯式Runge-Kutta法，簡稱R-K方法.它每步計(jì)算r個(gè)f值(即)，而ki由前面(i-1)個(gè)已算出的表示，故公式是顯式的.例如當(dāng)r=2時(shí)，公式可表示為 　

3、　　　　　　　　 　　　(3.2.5) 其中.改進(jìn)Euler法(3.1.11)就是一個(gè)二級顯式R-K方法.參數(shù)取不同的值，可得到不同公式. 3.2.2 二、三級顯式R-K方法 對r=2的顯式R-K方法(3.2.5)，要求選擇參數(shù)，使公式的精度階p盡量高，由局部截?cái)嗾`差定義　 (3.2.6) 令，對(3.2.6)式在處按Taylor公式展開，由于 將上述結(jié)果代入(3.2.6)得 要使公式(3.2.5)具有的階p=2,即，必須 　　　　　　　　　 　　(3.2.7) 即由此三式求的解不唯一.因r=2，由（3.2.5）式可知，于是有解 　　　　

4、　　　　　　　　　　 (3.2.8) 它表明使(3.2.5)具有二階的方法很多，只要都可得到二階精度R-K方法.若取，則，則得改進(jìn)Euler法(3.1.11)，若取，則得，此時(shí)(3.2.5)為 　　　　　　　　　　　　　　(3.2.9) 其中稱為中點(diǎn)公式.改進(jìn)的Euler法(3.1.11)及中點(diǎn)公式(3.2.9)是兩個(gè)常用的二級R-K方法，注意二級R-K方法只能達(dá)到二階，而不可能達(dá)到三階.因?yàn)閞=2只有4個(gè)參數(shù)，要達(dá)到p=3則在(3.2.6)的展開式中要增加3項(xiàng)，即增加三個(gè)方程,加上(3.2.7)的三個(gè)方程，共計(jì)六個(gè)方程求4個(gè)待定參數(shù)，驗(yàn)證得出是無解的.當(dāng)然r=2,p=

5、2的R-K方法(3.2.5)當(dāng)取其他數(shù)時(shí)，也可得到其他公式，但系數(shù)較復(fù)雜，一般不再給出. 對r=3的情形，要計(jì)算三個(gè)k值，即 其中　 將按二元函數(shù)在處按Taylor公式展開，然后代入局部截?cái)嗾`差表達(dá)式，可得 可得三階方法，其系數(shù)共有8個(gè)，所應(yīng)滿足的方程為 　　　　　　　　　 (3.2.10) 這是8個(gè)未知數(shù)6個(gè)方程的方程組，解也是不唯一的，通常.一種常見的三級三階R-K方法是下面的三級Kutta方法： 　　　　　　　　　　　 　　　(3.2.11) 附：R-K的三級Kutta方法程序如下 function y =

6、DELGKT3_kuta(f, h,a,b,y0,varvec) format long; N = (b-a)/h; y = zeros(N+1,1); y(1) = y0; x = a:h:b; var = findsym(f); for i=2:N+1 K1 = Funval(f,varvec,[x(i-1) y(i-1)]); K2 = Funval(f,varvec,[x(i-1)+h/2 y(i-1)+K1*h/2]); K3 = Funval(f,varvec,[x(i-1)+h y(i-1)-h*K1+K2*2*h]); y(i) =

7、 y(i-1)+h*(K1+4*K2+K3)/6; %滿足c1+c2+c3=1,(1/6 4/6 1/6) end format short; 3.2.3 四階R-K方法及步長的自動選擇 利用二元函數(shù)Taylor展開式可以確定(3.2.4)中r=4,p=4的R-K方法，其迭代公式為 其中，，而 共計(jì)13個(gè)參數(shù)待定，Taylor展開分析局部截?cái)嗾`差，使得精度達(dá)到四階，即誤差為。 于是，r=4，p=4的13個(gè)參數(shù)(c4不能為0)引出了多種方案和挑戰(zhàn)，如： 參數(shù)優(yōu)化使階數(shù)增加到5階，得到四階五階R-K方法，matlab中有程序ode45； 四級四階R-K方法的步長自動選

8、??； 結(jié)合新算法的應(yīng)用算法構(gòu)造； 適應(yīng)于新的領(lǐng)域?qū)崿F(xiàn)求解； …… 經(jīng)典的四階R-K方法是： 　　　　　　　　　 　　(3.2.12) 其中也需滿足，這里為(1/6 2/6 2/6 1/6).它的局部截?cái)嗾`差，故p=4，這是最常用的四階R-K方法，數(shù)學(xué)庫中都有用此方法求解初值問題的軟件.這種方法的優(yōu)點(diǎn)是精度較高，缺點(diǎn)是每步要算4個(gè)右端函數(shù)值，計(jì)算量較大. 例3.3 　用經(jīng)典四階R-K方法解例1.1的初值問題，仍取h=0.1，計(jì)算到，并與改進(jìn)Euler法、梯形法在處比較其誤差大小. 解 用四階R-K方法公式(3.2.12)，此處，于是當(dāng)n=0時(shí) 于是，按公式(3

9、.2.12)可算出 此方法誤差：　　　　　　　　　 改進(jìn)Euler法誤差：　　　　　　　 梯形法誤差：　　　　　　　　　 可見四階R-K方法的精度比二階方法高得多.用四階R-K方法求解初值問題(1.2.1)精度較高，但要從理論上給出誤差的估計(jì)式則比較困難.那么應(yīng)如何判斷計(jì)算結(jié)果的精度以及如何選擇合適的步長h？通常是通過不同步長在計(jì)算機(jī)上的計(jì)算結(jié)果近似估計(jì).設(shè)在處的值,當(dāng)時(shí)，的近似為，于是由四階R-K方法有 若以為步長，計(jì)算兩步到，則有 　　　　　　　　　 于是得 即 或 　　　　　　　　　　 　(3.2.13) 它給出了

10、誤差的近似估計(jì).如果(ε為給定精度)，則認(rèn)為以為步長的計(jì)算結(jié)果滿足精度要求，若，則還可放大步長.因此(3.2.13)提供了自動選擇步長的方法. 附：經(jīng)典的4級(也是4階)R-K方法的程序： function y = DELGKT4_Rungkuta(f, h,a,b,y0,varvec) format long; N = (b-a)/h; y = zeros(N+1,1); y(1) = y0; x = a:h:b; var = findsym(f); for i=2:N+1 K1 = Funval(f,varvec,[x(i-1) y(i-1)]); K

11、2 = Funval(f,varvec,[x(i-1)+h/2 y(i-1)+K1*h/2]); K3 = Funval(f,varvec,[x(i-1)+h/2 y(i-1)+K2*h/2]); K4 = Funval(f,varvec,[x(i-1)+h y(i-1)+h*K3]); y(i) = y(i-1)+h*(K1+2*K2+2*K3+K4)/6; end format short; 講解： 求初值問題（1.2.1）的單步法主要是指Runge-Kutta法，本節(jié)主要討論顯式R－K方法，建立具體的計(jì)算公式使用的是Taylor展開，

12、形如（3.2.4）的顯式R－K方法，當(dāng)r＝1時(shí)就是Euler法，因此只要討論的計(jì)算公式，在r確定后如何推導(dǎo)公式都是一樣的，只是r越大計(jì)算越復(fù)雜，為了掌握了解公式來源，只要以r＝2為例推導(dǎo)計(jì)算公式即可.因此本節(jié)重點(diǎn)就是用Taylor展開求出r＝2的顯式R-K方法的計(jì)算公式，由于方法的局部截?cái)嗾`差為（3.2.6），的右端有的項(xiàng)，要對它做Taylor展開，就要用到二元函數(shù)的Taylor展開，按照二元函數(shù)Taylor級數(shù) 　?。?.2.14） 將它用到（3.2.6）的的展開式中，即可得到按升冪整理出的結(jié)果，對r＝2的公式只能得到＝2階的公式，即，于是2級R-K方法（3.2.5）的系數(shù)必須滿足（3.

13、2.7）給出的方程，它的解由（3.2.8）給出，只要，求出的公式都是r=2的2階R-K方法.而常用的就是得到的改進(jìn)Euler法（3.1.11）和得到的中點(diǎn)公式（3.2.9）. 我們知道，顯式單步法的數(shù)值公式為： 如果取定右邊的 則稱該數(shù)值公式為顯式Runge-Kutta方法． （一）顯式Runge-Kutta方法 1. 二階二級R-K方法　 ， 其中，．四參數(shù)滿足 舉例：中點(diǎn)公式　　四個(gè)參數(shù)，即 ， 程序如下 function y = DELGKT2_mid(f, h,a,b,y0,varvec) format long; N = (b-a)/h;

14、y = zeros(N+1,1); y(1) = y0; x = a:h:b; var = findsym(f); for i=2:N+1 v1 = Funval(f,varvec,[x(i-1) y(i-1)]); t = y(i-1) + h*v1/2; v2 = Funval(f,varvec,[x(i-1)+h/2 t]); y(i) = y(i-1)+h*v2; end format short; 休恩（suen）法 c1=1/4, c2=3/4, a2=b21=2/3 function y = DELGKT2_suen(f,

15、h,a,b,y0,varvec) format long; N = uint16((b-a)/h); y = zeros(N+1,1); y(1) = y0; x = a:h:b; var = findsym(f); for i=2:N+1 K1 = Funval(f,varvec,[x(i-1) y(i-1)]); K2 = Funval(f,varvec,[x(i-1)+2*h/3 y(i-1)+K1*2*h/3]); y(i) = y(i-1)+h*(K1+3*K2)/4; end format short; 改進(jìn)的Euler法 c1

16、=1/2, c2=1/2, a2=b21=1 function y = DEModifEuler(f, h,a,b,y0,varvec) format long; N = (b-a)/h; y = zeros(N+1,1); y(1) = y0; x = a:h:b; var = findsym(f); for i=2:N+1 v1 = Funval(f,varvec,[x(i-1) y(i-1)]); t = y(i-1) + h*v1; v2 = Funval(f,varvec,[x(i) t]); y

17、(i) = y(i-1)+h*(v1+v2)/2; end format short; 實(shí)驗(yàn)對比參數(shù)取法不同得到不同的二階R-K方法精確程度各異：這個(gè)已經(jīng)可以做(同學(xué)們已做？！) 注意：二級R-K方法只能達(dá)到二階，而不可能達(dá)到三階.因?yàn)閞=2只有4個(gè)參數(shù)，要達(dá)到p=3則在(3.2.6)的展開式中要增加3項(xiàng)，即增加三個(gè)方程,加上(3.2.7)的三個(gè)方程，共計(jì)六個(gè)方程求4個(gè)待定參數(shù)，驗(yàn)證得出是無解的. 2．三階三級R-K方法 其中，，而 共計(jì)8個(gè)參數(shù)待設(shè)定。Taylor展開分析局部截?cái)嗾`差，使得精度達(dá)到三級，即誤差為。 于是，r=3，p=3的8個(gè)參數(shù)所滿足的關(guān)系

18、式(c3不能為0)為 Kutta法 function y = DELGKT3_kutta(f, h,a,b,y0,varvec) format long; N = (b-a)/h; y = zeros(N+1,1); y(1) = y0; x = a:h:b; var = findsym(f); for i=2:N+1 K1 = Funval(f,varvec,[x(i-1) y(i-1)]); K2 = Funval(f,varvec,[x(i-1)+h/2 y(i-1)+K1*h/2]); K3 = Funval(f,varvec,[x(i-1

19、)+h y(i-1)-h*K1+K2*2*h]); y(i) = y(i-1)+h*(K1+4*K2+K3)/6; end format short; Suen法 function y = DELGKT3_suen(f, h,a,b,y0,varvec) format long; N = (b-a)/h; y = zeros(N+1,1); y(1) = y0; x = a:h:b; var = findsym(f); for i=2:N+1 K1 = Funval(f,varvec,[x(i-1) y(i-1)]); K2 = Funva

20、l(f,varvec,[x(i-1)+h/3 y(i-1)+K1*h/3]); K3 = Funval(f,varvec,[x(i-1)+2*h/3 y(i-1)+K2*2*h/3]); y(i) = y(i-1)+h*(K1+3*K3)/4; end format short; 3. 四階四級R-K方法 舉例：經(jīng)典四級Rung-Kutta法 程序 function y = DELGKT4_Rungkutta(f, h,a,b,y0,varvec) format long; N = (b-a)/h; y = zeros(N+1,1); y(1

21、) = y0; x = a:h:b; var = findsym(f); for i=2:N+1 K1 = Funval(f,varvec,[x(i-1) y(i-1)]); K2 = Funval(f,varvec,[x(i-1)+h/2 y(i-1)+K1*h/2]); K3 = Funval(f,varvec,[x(i-1)+h/2 y(i-1)+K2*h/2]); K4 = Funval(f,varvec,[x(i-1)+h y(i-1)+h*K3]); y(i) = y(i-1)+h*(K1+2*K2+2*K3+K4)/6;

22、end format short; 基爾(jer、Gear)法 function y = DELGKT4_jer(f, h,a,b,y0,varvec) format long; N = (b-a)/h; C2 = sqrt(2); y = zeros(N+1,1); y(1) = y0; x = a:h:b; var = findsym(f); for i=2:N+1 K1 = Funval(f,varvec,[x(i-1) y(i-1)]); K2 = Funval(f,varvec,[x(i-1)+h/2 y(i-1)+K1*h/2]); K3 = F

23、unval(f,varvec,[x(i-1)+h/2 y(i-1)+K1*h*(C2-1)/2+K2*h*(2-C2)/2]); K4 = Funval(f,varvec,[x(i-1)+h y(i-1)-K2*h*C2/2+K3*h*(2+C2)/2]); y(i) = y(i-1)+h*(K1+(2-C2)*K2+(2+C2)*K3+K4)/6; end format short; 變型Rung-Kutta法 function y = DELGKT4_qt(f, h,a,b,y0,varvec) format long; N = (b-a)/h; y = zero

24、s(N+1,1); y(1) = y0; x = a:h:b; var = findsym(f); for i=2:N+1 K1 = Funval(f,varvec,[x(i-1) y(i-1)]); K2 = Funval(f,varvec,[x(i-1)+h/3 y(i-1)+K1*h/3]); K3 = Funval(f,varvec,[x(i-1)+2*h/3 y(i-1)-K1*h/3+K2*h]); K4 = Funval(f,varvec,[x(i-1)+h y(i-1)+h*(K1-K2+K3)]); y(i) = y(i-1)+

25、h*(K1+3*K2+3*K3+K4)/8; end format short; 說明： 四階四級R-K方法的步長自動選取問題 用四階R-K方法求解初值問題(1.2.1)精度已經(jīng)足夠高了，但要從理論上給出誤差的估計(jì)式則比較困難.那么應(yīng)如何判斷計(jì)算結(jié)果的精度以及如何選擇合適的步長h？通常是通過不同步長在計(jì)算機(jī)上的計(jì)算結(jié)果近似估計(jì).設(shè)在處的值,當(dāng)時(shí)，的近似為，于是由四階R-K方法有 若以為步長，計(jì)算兩步到，則有 　　　　　　　　　 于是得 即 或 　　　　　　　　　　 　(3.2.13) 它給出了誤差的近似估計(jì).如果 (ε

26、為給定精度，如1E-3)，則認(rèn)為以為步長的計(jì)算結(jié)果滿足精度要求，若，則還可放大步長.因此(3.2.13)提供了自動選擇步長的方法. （二）隱式(半隱式)R-K方法舉例 一種半隱式舉例，該格式可以通過除法化為顯式單步法格式，受到除法復(fù)雜度限制，容易造成耗時(shí)，不穩(wěn)，相對精度降低。 羅賽布諾克半隱式法 function y = DELSBRK(f, h,a,b,y0,varvec) format long; a1 = 1.40824829;a2 = 0.59175171; c1 = 0.17378667;c2 = c1; w1 = -0.41315432;w2 = 1.413

27、15432; N = (b-a)/h; y = zeros(N+1,1); y(1) = y0; x = a:h:b; var = findsym(f); dy = diff(f, varvec(2)); for i=2:N+1 f1 = Funval(f,varvec,[x(i-1) y(i-1)]); dy1 = Funval(dy,varvec,[x(i-1) y(i-1)]); k1 = h*f1/(1-h*a1*dy1); dy2 = Funval(dy,varvec,[x(i-1)+c1*h y(i-1)+c2*k1]);

28、 f2 = Funval(f,varvec,[x(i-1)+c1*h y(i-1)+c2*k1]); k2 = h*f2/(1-h*a2*dy2); y(i) = y(i-1)+w1*k1+w2*k2; end format short; 3.3 單步法的收斂性與絕對穩(wěn)定性 3.3.1 單步法的收斂性 定義3.1 設(shè)y(x)是初值問題(1.2.1)的精確解，是單步法(3.2.2)在處產(chǎn)生的近似解，若 則稱方法(3.2.2)產(chǎn)生的數(shù)值解收斂于. 實(shí)際上，定義中是一固定點(diǎn)，當(dāng)h→0時(shí)n→∞，n不是固定的.因顯然方法收斂，則在固定點(diǎn)處的整體誤差，當(dāng)

29、p≥1時(shí). 下面定理給出方法(3.2.2)收斂的條件. 定理3.1　設(shè)初值問題(1.2.1)的單步法(3.2.2)是p階方法(p≥1)，且函數(shù)對y滿足Lipschitz條件，即存在常數(shù)L＞0，使對，均有 　　　 　 則方法(3.2.2)收斂，且. 定理證明略. 3.3.2 絕對穩(wěn)定性 用單步法(3.2.2)求數(shù)值解，由于原始數(shù)據(jù)及計(jì)算過程舍入誤差影響，實(shí)際得到的不是而是，其中是誤差，再計(jì)算下一步得到 以Euler法為例，若令,則 　　　　(3.3.1) 如果，則從計(jì)算到誤差不增長，它是穩(wěn)定的.但如果條件不滿足就不穩(wěn)定. 例3.4 y′=-100y,y(0)=1,精

30、確解為，用Euler法求解得 若取h=0.025，則,當(dāng)，而，顯然計(jì)算是不穩(wěn)定的. 如果用后退Euler法(3.1.5)解此例，仍取h=0.025,則 　　　　　　,即 顯然當(dāng)，計(jì)算是穩(wěn)定的. 　　由此看到穩(wěn)定性與方法有關(guān)，也與有關(guān)，在此例中.在研究方法的穩(wěn)定性時(shí)，通常不必對一般的f(x,y)進(jìn)行討論，而只針對模型方程 　　　　　　　　　　　　 　　(3.3.2) 這里可能為復(fù)數(shù).規(guī)定是因?yàn)闀r(shí)微分方程(3.3.2)本身是不穩(wěn)定的，而討論數(shù)值方法(3.2.2)的穩(wěn)定性，必須在微分方程本身穩(wěn)定的前提下進(jìn)行.另一方面，對初值問題(1.2.1)，若將f(x,y)在

31、處線性展開，可得 于是方程(1.2.1)可近似表示為 它表明用模型方程(3.3.2)是合理的，至于模型方程(3.3.2)中所以用復(fù)數(shù)λ是因?yàn)槌踔祮栴}(1.2.1)如果是方程組，即，則是(m×m)階矩陣，其特征值可能是復(fù)數(shù).當(dāng)然對單個(gè)方程，λ就是實(shí)數(shù)，此時(shí)只要規(guī)定＜0即可. 用單步法(3.2.2)解模型方程(3.3.2)可得到 　　　　　　　　　　　　　　　　 　(3.3.3) 其中依賴所選方法，如用Euler法，則 　　　　　　　　　　(3.3.4) 此時(shí)由(3.3.1)看到誤差方程也為，與(3.3.4)是一樣的.因此對一般單步法(3.2.2)誤差方程也

32、與(3.3.3)一致.下面再考慮二階R-K方法有 對四階R-K方法，可得 定義3.2　將單步法(3.2.2)用于解模型方程(3.3.2)，若得到(3.3.3)中的 則稱方法是絕對穩(wěn)定的.在復(fù)平面上復(fù)變量滿足 的區(qū)域，稱為方法(3.2.2)的絕對穩(wěn)定域，它與實(shí)軸的交點(diǎn)稱為絕對穩(wěn)定區(qū)間. 例如對Euler法， 在復(fù)平面上是以(-1，0)為圓心，以1為半徑的單位圓域內(nèi)部，當(dāng)為實(shí)數(shù)時(shí)，則得絕對穩(wěn)定區(qū)間為，因＜0，故有.在例3.4中 時(shí)方法穩(wěn)定，而例中取h=0.025故不穩(wěn)定. 對后退Euler法(3.1.5)， 因＜0，故，其絕對穩(wěn)定域是以(1，0)為圓心的單位圓外部，絕對穩(wěn)定區(qū)

33、間為，即對任何h＞0方法都是絕對穩(wěn)定的. 二階R-K方法的絕對穩(wěn)定區(qū)間為. 三階R-K方法的絕對穩(wěn)定區(qū)間為. 四階R-K方法的絕對穩(wěn)定區(qū)間為. 例3.5 用經(jīng)典四階R-K方法計(jì)算初值問題 　　　　　　　　　　 步長取h=0.1及0.2,給出計(jì)算誤差并分析其穩(wěn)定性. 解 本題直接按R-K方法(3.2.12)的公式計(jì)算.因精確解為，其計(jì)算誤差如表所示. 從計(jì)算結(jié)果看到，h=0.2時(shí)誤差很大，這是由于在λ=-20,h=0.2時(shí)λh=-4，而四階R-K方法的絕對穩(wěn)定區(qū)間為［-2.785,0］,故h=0.2時(shí)計(jì)算不穩(wěn)定，誤差很大.而h=0.1時(shí)=-2，其值在絕對穩(wěn)定區(qū)間［-2.78

34、5,0］內(nèi)，計(jì)算穩(wěn)定，故結(jié)果是可靠的. 講解： 由于微分方程初值問題數(shù)值解公式求出的解是一個(gè)逐次遞推的過程，因此原始數(shù)據(jù)誤差及計(jì)算過程舍入誤差對解的影響就是數(shù)值方法絕對穩(wěn)定性研究的問題，如果由計(jì)算誤差不增長，方法就是絕對穩(wěn)定的.為使問題得到簡化通常就是將方法用于解模型方程（3.3.2），對于單步法得到的差分方程為，由于模型方程的，代入Euler法，得，對二階R-K方法，例如，用改進(jìn)Euler法 于是　　　　　 對三階R-K方法有 對四階R-K方法有 只要方法，就是絕對穩(wěn)定的，這時(shí)的值當(dāng)n增大式是減少的，故計(jì)算穩(wěn)定.這時(shí)舍入誤差影響可忽略不計(jì)，而當(dāng)，則增大，方法不穩(wěn)定，計(jì)算

35、結(jié)果是不可靠的.因此用顯式單步法必須使，也就是步長選擇要滿足這一要求. 對于隱式的梯形公式 將模型方程，即代入得 于是　　注意，于是有 　　　　，對成立. 這就表明對任意步長h，梯形法都是絕對穩(wěn)定的. 3.4 練習(xí)題 2. 對于一階微分方程初值問題，用Euler法，改進(jìn)Euler法，二階R-K方法求解，并作圖比較。 3. 對于一階微分方程初值問題，用Euler法，改進(jìn)Euler法，二階R-K法，三階R-K法求解，并與真解作圖比較，列出誤差對比表格。 4. 對于一階微分方程初值問題，用三階、四階R-K法與真解作圖對比。 5. 對于一階微分方程初值問題，R-K法的半隱式、四階與真解作圖對比。 6. 對一階微分方程初值問題進(jìn)行穩(wěn)定性分析： a. 得出二階R-K法中步長的穩(wěn)定區(qū)間(域)； b. 自選兩個(gè)步長(h1為區(qū)間內(nèi)的數(shù)(如0.1), h2為區(qū)間外的數(shù)(如0.3))得到兩個(gè)不同數(shù)值解和精確解作圖比較，并列出誤差表格。