微分方程常用的兩種數(shù)值解法：歐拉方法和龍格庫塔法

1、 PAGE1 / NUMPAGES33 某師X大學(xué)本科畢業(yè)論文微分方程常用的兩種數(shù)值解法：歐拉方法與龍格庫塔法學(xué)生某XXX院系名稱數(shù)學(xué)與軟件科學(xué)學(xué)院專業(yè)名稱信息與計算科學(xué)班 級2006級 4 班學(xué) 號20060640XX指導(dǎo)教師Xxx某師X大學(xué)教務(wù)處二一年五月微分方程常用的兩種數(shù)值解法：歐拉方法與龍格庫塔法學(xué)生某：xxx 指導(dǎo)教師：xx【內(nèi)容摘要】微分方程是最有生命力的數(shù)學(xué)分支，在自然科學(xué)的許多領(lǐng)域中，都會遇到常微分方程的求解問題。當(dāng)前計算機的發(fā)展為常微分方程的應(yīng)用及理論研究提供了非常有力的工具，利用計算機解微分方程主要使用數(shù)值方法，歐拉方法和龍格庫塔方法是求解微分方程最典型常用的數(shù)值方法。本

2、文詳細研究了這兩類數(shù)值計算方法的構(gòu)造過程，分析了它們的優(yōu)缺點，以及它們的收斂性，相容性，及穩(wěn)定性。討論了步長的變化對數(shù)值方法的影響和系數(shù)不同的同階龍格庫塔方法的差別。通過編制C程序在計算機上實現(xiàn)這兩類方法及對一些典型算例的結(jié)果分析比較，能更深切體會它們的功能，優(yōu)缺點及適用場合，從而在實際應(yīng)用中能對不同類型和不同要求的常微分方程會選取適當(dāng)?shù)那蠼夥椒āｊP(guān)鍵詞：顯式單步法 歐拉（Euler）方法 龍格庫塔（RungeKutta）方法 截斷誤差 收斂性Two monly used numerical solution of differential equations：Euler method and

3、 Runge - Kutta methodStudent Name: Xiong Shiying Tutor：Zhang Li【Abstract】The differential equation is the most vitality branch in mathematics. In many domains of natural science, we can meet the ordinary differentialequation solution question. Currently, the development of puter has provided the ext

4、remely powerful tool for the ordinary differential equation application and the fundamental research, the puter solving differential equation mainly uses value method.The Euler method and the RungeKutta method are the most typical monly value method to solve the differential equation. This article d

5、issects the structure process of these two kinds of values monly value method to solve the analyses their good and bad points, to their astringency, the patibility, and the stability has made the proof. At the same time, the article discuss the length of stride to the numerical method changing influ

6、ence and the difference of the coefficient different same step Rungekutta method. Through establishing C program on the puter can realize these two kind of methods, Anglicizing some models of calculate example result can sincerely realize their function, the advantage and disadvantage points and the

7、 suitable situation, thus the suitable solution method can be selected to solve the different type and the different request ordinary differential equation in the practical application .Keywords:Explicit single-step processEuler methodRungeKutta methodtruncation errorconvergence目 錄微分方程常用的兩種數(shù)值解法：歐拉方法

8、與龍格庫塔法前言 常微分方程的形成與發(fā)展是和力學(xué)、天文學(xué)、物理學(xué)以及其他科學(xué)技術(shù)的發(fā)展密切相關(guān)的。數(shù)學(xué)其他分支的新發(fā)展，如復(fù)變函數(shù)、群、組合拓?fù)鋵W(xué)等，都對常微分方程的發(fā)展產(chǎn)生了深刻的影響，當(dāng)前計算機的發(fā)展更是為常微分方程的應(yīng)用及理論研究提供了非常有力的工具。牛頓在研究天體力學(xué)和機械力學(xué)的時候，利用了微分方程這個工具，從理論上得到了行星運動規(guī)律。后來，法國天文學(xué)家勒維烈和英國天文學(xué)家亞當(dāng)斯使用微分方程各自計算出那時尚未發(fā)現(xiàn)的海王星的位置。這些都使數(shù)學(xué)家更加深信微分方程在認(rèn)識自然、改造自然方面的巨大力量，微分方程也就成了最有生命力的數(shù)學(xué)分支。然而，我們知道，只有少數(shù)十分簡單的微分方程能夠用初等方法

9、求得它們的解，多數(shù)情形只能利用近似方法求解。在常微分方程課中的級數(shù)解法，逐步逼近法等就是近似解法。這些方法可以給出解的近似表達式，通常稱為近似解析方法。還有一類近似方法稱為數(shù)值方法，它可以給出解在一些離散點上的近似值，利用計算機解微分方程主要使用數(shù)值方法。本文主要討論一階常微分方程初值問題 （1.1）在區(qū)間上的數(shù)值解法，其中為關(guān)于，的已知函數(shù)，為給定的初始值，將上述問題的精確解記為。該問題常用的數(shù)值解法有：歐拉（Euler）方法、龍格庫塔（RungeKutta）方法及一些常用的線性多步法。本文重點介紹歐拉（Euler）方法和龍格庫塔（RungeKutta）方法。并對這兩種方法編制程序，體會它們

10、的功能、優(yōu)缺點及適用場合，對不同類型常微分方程會選取適當(dāng)?shù)那蠼夥椒ā?基本概念和準(zhǔn)備知識一階常微分方程初值問題是：其中是平面上某個區(qū)域上的連續(xù)函數(shù)，式（1.1.1）的微分方程一般有無窮多個解，式（1.1.2）是確定解的初始條件，如果一元函數(shù)對一切滿足（1）；（2）；（3）存在而且；則稱是初值問題（1.1）在區(qū)間上的解。誤差：假定在計算時，用到的前一步的值是準(zhǔn)確的即，把用計算得到的近似值記為，估計誤差：= y(xn+1)，這種誤差稱為局部截斷誤差。如果不作這一假定，在每一步計算時除局部截斷誤差以外，還有由于前一步不準(zhǔn)確而引起的誤差，稱為總體截斷誤差。收斂性：對于解初值問題的數(shù)值方法，我們希望它產(chǎn)

11、生的數(shù)值解收斂于初值問題的準(zhǔn)確解，“收斂性”的一般定義為：對于所有滿足引理1.1條件的初值問題（1.1），如果有一種顯式單步法：產(chǎn)生的近似解，對于任意固定的，均有，則稱該顯式單步法是收斂的。相容性：顯式單步法（1.2.1）稱為與原微分方程相容，如果 （1.2.3）成立。并稱式（1.2.3）為相容性條件。穩(wěn)定性：在實際計算中，一方面初始值不一定精確，往往帶有一定誤差，同時由于計算機字長有限，在計算過程中有舍入誤差，而且這種誤差在式（1.2.1）的遞推過程中會傳遞下去，對以后的結(jié)果產(chǎn)生影響。因此要考慮舍入誤差的積累是否會得到控制，也即要考慮數(shù)值方法的穩(wěn)定性。當(dāng)時，若舍入誤差引起的后果是有限的，則可

12、以認(rèn)為該方法是數(shù)值穩(wěn)定的。2 歐拉方法2.1歐拉方法簡介對常微分方程初值問題（1.1）用數(shù)值方法求解時，我們總是認(rèn)為（1.1）的解存在且唯一。歐拉方法是解初值問題的最簡單、最原始的數(shù)值方法，它是顯式單步法。下面介紹幾種導(dǎo)出歐拉法的途徑，每個途徑皆可以推導(dǎo)出更為有效的數(shù)值方法。（1）Taylor展開在點將作Taylor展開得：當(dāng)充分小時略去誤差項，并注意到，得，以近似代替，以近似代替，且用“=”代替“”得差分方程初值問題：， (2.1.1)用式(2.1.1)求解初值問題（1.1）的方法稱為歐拉方法。（2）數(shù)值微分由導(dǎo)數(shù)的定義知，對于充分小的整理得，對此作相似的處理也可以得到歐拉方法(2.1.1)

13、。（3）數(shù)值積分在區(qū)間對積分得 （2.1.2）用數(shù)值積分的左矩形公式計算式（2.1.2）右端的積分，得，于是同樣可以得到歐拉方法(2.1.1)。（4）多項式插值利用解和其導(dǎo)數(shù)在點的值，作一次埃爾米特插值，得到關(guān)于的插值多項式：，用近似代替就得到歐拉方法。（5）待定系數(shù)法在第步，已知和，利用這兩個值估計出下一步的，將已知的值與估計值作線性組合：，其中， 為待定系數(shù)。為確定這兩個參數(shù)，要求這個估計值對和（為常數(shù)）精確成立。如果，則，得到方程：，得。如果，則，這樣有：，說明，這樣估計值為：，即為歐拉方法。歐拉方法的幾何意義：由點斜式得切線方程等步長為，則，可由切線方程算出：，逐步計算出在點的值：，用

14、分段的折線逼近函數(shù)為“折線法”而非“切線法”，除第一個點是曲線上的切線，其它都不是。0圖1歐拉方法的幾何意義2.2歐拉方法的截斷誤差，收斂性，相容性，穩(wěn)定性設(shè)，把在處展開成泰勒級數(shù)，即再由歐拉方法： 兩式相減得歐拉方法的局部截斷誤差為：，若在上充分光滑，且令，則，故歐拉方法是一階方法或具有一階精度。歐拉方法的增量函數(shù)就是，由引理1.3、引理1.4知當(dāng)滿足Lipschitz條件時歐拉方法是收斂的而且是相容的。用歐拉方法求解典型方程（1.2.4）的計算公式為：，有 。要讓，必須有，因此歐拉方法的絕對穩(wěn)定域為 ：，當(dāng)為實數(shù)時，絕對穩(wěn)定區(qū)間為。在復(fù)平面上，是以1為半徑、以為圓心的內(nèi)部。3 龍格庫塔法3

15、.1 龍格庫塔法的基本思想為了導(dǎo)出龍格庫塔法的一般公式，我們?nèi)∪缦碌木€性組合形式： （3.2.3）其中 （3.2.4）即 ,；a21,a31除c1=0外均為待定系數(shù)。顯然用公式（3.2.3）每計算一個新值要計算函數(shù)的值s次，又因每個都能以一種明顯的方式由,計算出來，故將公式（3.2.3）稱為s級顯式龍格庫塔法。s級顯式龍格庫塔法又可以寫成下面既簡潔又直觀的陣列形式： 03.2 二、三、四級龍格庫塔法 常見的二級二階龍格庫塔法有：中點方法（?。┒AHeun方法（取）取s = 3，完全仿照上述方法推導(dǎo)出三階龍格庫塔公式。這時參數(shù)滿足下列條件這個方程組解也不唯一，從而可以得到不同的三級三階龍格庫塔法

16、。兩個較為常見的三級三階龍格庫塔法是：Kutta方法（取=， = 1，a21 = ，a31 = ，a32 = 2， = ， = ， = ）將它代入（3.2.3）得(3.3.2)三階Heun方法（取=， = ，a21 = ，a31 = 0，a32 = ， = ， = ， =）將它代入（3.2.3）得通常人們所說的龍格庫塔法是指四階而言的。取s=4，我們同樣可以仿照二階的情形推導(dǎo)出此公式。這樣就得到如下含12個未知數(shù)但僅有10個方程的方程組此方程組的解同樣不是唯一的，從而有許多不同的四級四階龍格庫塔法，最常用的兩個四階公式是：經(jīng)典龍格庫塔法 （3.3.3）RKG(Runge-Kutta-Gill)

17、方法 (3.3.4)經(jīng)典龍格庫塔方法的幾何解釋：將（3.3.3）式中的改記為，從而得到經(jīng)典龍格庫塔方法的計算公式為只考慮區(qū)間上的解曲線。由(3.3.3)式可知，曲線在處的斜率，在處的斜率，而和則是在中點處的斜率的兩個近似值，如圖2：圖2另一方面有：，用Simpson積分公式逼近此式右端的積分，得上式右端要計算的三個值。由前面的結(jié)果，可取，見圖3：圖3代入后得到 （3.3.5）另外，從(3.3.3)式容易看出，當(dāng)與無關(guān)時，公式(3.3.5)實際上就是用標(biāo)準(zhǔn)的Simpson積分公式計算積分得到的，因此，可以將經(jīng)典龍格庫塔方法視為是用Simpson積分公式的推廣形式計算積分而得到的數(shù)值方法。經(jīng)典龍格

18、庫塔法的穩(wěn)定性：增量函數(shù)：，其中代入后得于是有=從而得出龍格庫塔方法的絕對穩(wěn)定區(qū)間為 1經(jīng)典龍格庫塔方法的算法框圖：開 始讀入輸出=結(jié) 束圖4 經(jīng)典龍格庫塔法算法框圖4 應(yīng)用舉例建模及其分析4.1 歐拉方法解題及其數(shù)學(xué)模型4.1.1 問題提出例4.1.1假定某公司的凈資產(chǎn)因資產(chǎn)本身產(chǎn)生了利息而以4%的年利率增長，同時，該公司以每年100萬的數(shù)額支付職工工資。凈資產(chǎn)的微分方程為 （t以年為單位）分別以初始值問題：用Euler公式預(yù)測公司24后的凈資產(chǎn)趨勢。4.1.2模型建立 分析：這是求微分方程的數(shù)值積分，為的是預(yù)測公司24年后的凈資產(chǎn)趨勢。確定變量：設(shè)凈資產(chǎn)是時間的微分函數(shù)，不妨設(shè)變量t為時間

19、（以年為單位）。設(shè)為第n年的凈資產(chǎn)，為第n+1年的凈資產(chǎn)，利息以每年4%的速度增長，且公司每年支付職工工資為100萬，則第n+1年的凈資產(chǎn)增長數(shù)額為，由于已得第n年的凈資產(chǎn)為，則第n年的凈資產(chǎn)加上第n+1年的凈資產(chǎn)增長數(shù)額就得到第n+1年的凈資產(chǎn)。歸納公式： （4-1）確定其微分方程的標(biāo)準(zhǔn)形式，這就是例4.1.1的微分方程模型。4.1.3 解決問題分別以代入，x表示利息贏利低于工資支出的數(shù)額，y表示利息贏利與工資支出平衡的數(shù)額，z表示利息贏利高于工資支出的數(shù)額，計算結(jié)果見表1.表1nn11460.25003540.13834.92625004165.0721418.425003581.6147

20、68.32425004231.6831375.1425003624.8615699.05625004300.9441330.1425003669.8616627.01925004372.9851283.3525003716.6517552.125004447.961234.6825003765.32184748271184.0725003815.93193938581131.4325003868.5720308.87725004691.1291076.6925003923.3121221.23225004778.77101019.762500

21、3980.2422130.08125004869.9211960.54625004039.452335.284525004964.7212898.96825004101.0324-63.304225005063.3從表1可以看到當(dāng)利息贏利低于工資的支出，公司的凈資產(chǎn)逐年減少，以致凈資產(chǎn)為負(fù)值；當(dāng)利息贏利與工資的支出平衡時，公司的凈資產(chǎn)每年保持不變；當(dāng)利息贏利超過工資的支出，公司的凈資產(chǎn)穩(wěn)步增長。再用歐拉方法求解，每計算一步，依次需要計算1次、2次和4次函數(shù)的值，為了比較在計算量相同的條件下近似解的精度，步長分別取0.025、0.05、0.1。我對歐拉的C程序做了一些優(yōu)化，加入了計算誤差，計算結(jié)

22、果如下（表示近似解、表示準(zhǔn)確解、表示誤差）：This is Euler method:please input a b and y0：0 1 1 please input N：40 x0=0.000000 y0=1.000000 y10=1.000000 e0=0.000000 x1=0.025000 y1=1.025000 y11=1.024695 e1=-0.000305x2=0.050000 y2=1.049405 y12=1.048809 e2=-0.000597x3=0.075000 y3=1.073258 y13=1.072381 e3=-0.000878x4=0.100000 y

23、4=1.096596 y14=1.095445 e4=-0.001151x5=0.125000 y5=1.119451 y15=1.118034 e5=-0.001417x6=0.150000 y6=1.141854 y16=1.140175 e6=-0.001679x7=0.175000 y7=1.163832 y17=1.161895 e7=-0.001937x8=0.200000 y8=1.185410 y18=1.183216 e8=-0.002194x9=0.225000 y9=1.206609 y19=1.204159 e9=-0.002450 x10=0.250000 y10=1

24、.227451 y110=1.224745 e10=-0.002706x11=0.275000 y11=1.247954 y111=1.244990 e11=-0.002964x12=0.300000 y12=1.268134 y112=1.264911 e12=-0.003223x13=0.325000 y13=1.288009 y113=1.284523 e13=-0.003486x14=0.350000 y14=1.307593 y114=1.303841 e14=-0.003753x15=0.375000 y15=1.326900 y115=1.322876 e15=-0.004024

25、x16=0.400000 y16=1.345941 y116=1.341641 e16=-0.004300 x17=0.425000 y17=1.364730 y117=1.360147 e17=-0.004583x18=0.450000 y18=1.383278 y118=1.378405 e18=-0.004873x19=0.475000 y19=1.401594 y119=1.396424 e19=-0.005170 x20=0.500000 y20=1.419689 y120=1.414214 e20=-0.005475x21=0.525000 y21=1.437572 y121=1.

26、431782 e21=-0.005790 x22=0.550000 y22=1.455251 y122=1.449138 e22=-0.006113x23=0.575000 y23=1.472735 y123=1.466288 e23=-0.006447x24=0.600000 y24=1.490032 y124=1.483240 e24=-0.006792x25=0.625000 y25=1.507149 y125=1.500000 e25=-0.007149x26=0.650000 y26=1.524093 y126=1.516575 e26=-0.007518x27=0.675000 y

27、27=1.540872 y127=1.532971 e27=-0.007900 x28=0.700000 y28=1.557490 y128=1.549193 e28=-0.008297x29=0.725000 y29=1.573955 y129=1.565248 e29=-0.008708x30=0.750000 y30=1.590273 y130=1.581139 e30=-0.009134x31=0.775000 y31=1.606449 y131=1.596872 e31=-0.009577x32=0.800000 y32=1.622489 y132=1.612452 e32=-0.0

28、10037x33=0.825000 y33=1.638397 y133=1.627882 e33=-0.010515x34=0.850000 y34=1.654180 y134=1.643168 e34=-0.011013x35=0.875000 y35=1.669842 y135=1.658312 e35=-0.011530 x36=0.900000 y36=1.685388 y136=1.673320 e36=-0.012068x37=0.925000 y37=1.700823 y137=1.688194 e37=-0.012629x38=0.950000 y38=1.716151 y13

29、8=1.702939 e38=-0.013212x39=0.975000 y39=1.731376 y139=1.717556 e39=-0.013820 x40=1.000000 y40=1.746504 y140=1.732051 e40=-0.014453this is improve Euler method please input a b and y0:：0 1 1 please input N：20 x0=0.000000 y0=1.000000 y10=1.000000 e0=0.000000 x1=0.050000 y1=1.048869 y11=1.048809 e1=-0

30、.000060 x2=0.100000 y2=1.095561 y12=1.095445 e2=-0.000116x3=0.150000 y3=1.140345 y13=1.140175 e3=-0.000169x4=0.200000 y4=1.183437 y14=1.183216 e4=-0.000221x5=0.250000 y5=1.225017 y15=1.224745 e5=-0.000272x6=0.300000 y6=1.265236 y16=1.264911 e6=-0.000325x7=0.350000 y7=1.304219 y17=1.303841 e7=-0.0003

31、78x8=0.400000 y8=1.342074 y18=1.341641 e8=-0.000433x9=0.450000 y9=1.378896 y19=1.378405 e9=-0.000492x10=0.500000 y10=1.414766 y110=1.414214 e10=-0.000553x11=0.550000 y11=1.449756 y111=1.449138 e11=-0.000618x12=0.600000 y12=1.483927 y112=1.483240 e12=-0.000687x13=0.650000 y13=1.517337 y113=1.516575 e

32、13=-0.000762x14=0.700000 y14=1.550035 y114=1.549193 e14=-0.000841x15=0.750000 y15=1.582066 y115=1.581139 e15=-0.000927x16=0.800000 y16=1.613472 y116=1.612452 e16=-0.001020 x17=0.850000 y17=1.644289 y117=1.643168 e17=-0.001121x18=0.900000 y18=1.674551 y118=1.673320 e18=-0.001230 x19=0.950000 y19=1.70

33、4288 y119=1.702939 e19=-0.001349x20=1.000000 y20=1.733529 y120=1.732051 e20=-0.0014784.1.4 結(jié)論這道題用普通的微分方程也能列式求解，關(guān)鍵在于如何預(yù)測若干年后的凈資產(chǎn)趨勢，用普通的微分方程就無法進行預(yù)測，且人工計算量相當(dāng)大，這里我使用歐拉方法可以計算出精度以及誤差，通過電腦運行程序就可以預(yù)測出若干年后的凈資產(chǎn)情況，歐拉方法的使用使得解題更加方便且精確。4.2 龍格庫塔解題及其數(shù)學(xué)模型4.2.1 問題提出 例4.2.1 兩種果樹寄生蟲，其數(shù)量分別是其中一種寄生蟲以吃另一種寄生蟲為生，兩種寄生蟲的增長函數(shù)如下列

34、常微分方程組所示：問題：預(yù)測3年后這一對寄生蟲的數(shù)量。4.2.2 模型建立分析：這是一個典型的常微分方程例題，要求預(yù)測出3年后這對寄生蟲的數(shù)量。確定變量：假定時間是寄生蟲數(shù)量的積分函數(shù)，不妨設(shè)變量t為時間（以年為單位）。由題知有兩種寄生蟲u和v，u寄生蟲年增長函數(shù)為，v寄生蟲年增長函數(shù)為，初始值：u寄生蟲為1.6，v寄生蟲為1.2，由于其中一種寄生蟲以吃另一種寄生蟲為生，我們可建立u和v的關(guān)聯(lián)函數(shù)f(u,v),g(u,v)。歸納后得到公式： （4-2）（4-2）即為例4.2.1所述問題的微分方程模型。4.2.3 解決問題在本題中用Euler預(yù)估校正公式取h=1，計算結(jié)果如表2.表2t/年u(t

35、)v(t)11.61.221.024571.2683430.6409121.336640.3912111.41077我把龍格庫塔的C程序進行編輯后得到結(jié)果:（ei表示誤差）please input a b and y0:：0 1 1 please input N： 10 x0=0.000000 y0=1.000000 y10=1.000000 e0=0.000000 x1=0.100000 y1=1.095446 y11=1.095445 e1=-0.000000 x2=0.200000 y2=1.183217 y12=1.183216 e2=-0.000001x3=0.300000 y3=1

36、.264912 y13=1.264911 e3=-0.000001x4=0.400000 y4=1.341642 y14=1.341641 e4=-0.000001x5=0.500000 y5=1.414215 y15=1.414214 e5=-0.000002x6=0.600000 y6=1.483242 y16=1.483240 e6=-0.000002x7=0.700000 y7=1.549196 y17=1.549193 e7=-0.000003x8=0.800000 y8=1.612455 y18=1.612452 e8=-0.000004x9=0.900000 y9=1.67332

37、4 y19=1.673320 e9=-0.000004x10=1.000000 y10=1.732056 y110=1.732051 e10=-0.0000054.2.3 結(jié)論從上面的結(jié)果可以分析，用每一種方法計算節(jié)點=0.1、0.20.9、1.0，上的值都需要計算4次的值，即它們的計算量基本相同，其結(jié)果是經(jīng)典的龍格庫塔方法的精度最好，龍格庫塔方法的推導(dǎo)是基于Taylor級數(shù)方法的，因而在使用高階龍格庫塔方法計算時可以很精確的推算出寄生蟲每一年的數(shù)量。4.3 分別使用二階、三階龍格庫塔法解初值問題對一些特殊的微分方程，使用歐拉方法和低階的龍格庫塔方法也能達到很高的精度，例如：微分方程的解析解是

38、一次函數(shù)，則用歐拉方法求得的數(shù)值解與準(zhǔn)確解相符，微分方程的解析解是二次函數(shù)，則用二階龍格庫塔方法求得的數(shù)值解與準(zhǔn)確解相符。微分方程的解析解是三次多項式，用三階龍格庫塔方法求得的數(shù)值解與準(zhǔn)確解相符。4.3.1 建立模型建立初值問題模型： (4.1)4.3.2 用不同的龍格庫塔法解（4.1）初值問題模型 （1）用二階龍格庫塔方法求解初值問題（4.1）結(jié)果如下：this is the second-order Runge-Kutta(Heun) method please input a b and y0： 0 1 1 please input N：10 x0=0.000000 y0=1.00000

39、0 y10=1.000000 e0=0.000000 x1=0.100000 y1=1.110000 y11=1.110000 e1=0.000000 x2=0.200000 y2=1.240000 y12=1.240000 e2=0.000000 x3=0.300000 y3=1.390000 y13=1.390000 e3=0.000000 x4=0.400000 y4=1.560000 y14=1.560000 e4=0.000000 x5=0.500000 y5=1.750000 y15=1.750000 e5=0.000000 x6=0.600000 y6=1.960000 y16=

40、1.960000 e6=0.000000 x7=0.700000 y7=2.190000 y17=2.190000 e7=0.000000 x8=0.800000 y8=2.440000 y18=2.440000 e8=0.000000 x9=0.900000 y9=2.710000 y19=2.710000 e9=0.000000 x10=1.000000 y10=3.000000 y110=3.000000 e10=0.000000this is improve Euler methodplease input a b and y0：0 1 1 please input N：10 x0=0

41、.000000 y0=1.000000 y10=1.000000 e0=0.000000 x1=0.100000 y1=1.110000 y11=1.110000 e1=0.000000 x2=0.200000 y2=1.240000 y12=1.240000 e2=0.000000 x3=0.300000 y3=1.390000 y13=1.390000 e3=0.000000 x4=0.400000 y4=1.560000 y14=1.560000 e4=0.000000 x5=0.500000 y5=1.750000 y15=1.750000 e5=0.000000 x6=0.60000

42、0 y6=1.960000 y16=1.960000 e6=0.000000 x7=0.700000 y7=2.190000 y17=2.190000 e7=0.000000 x8=0.800000 y8=2.440000 y18=2.440000 e8=0.000000 x9=0.900000 y9=2.710000 y19=2.710000 e9=0.000000 x10=1.000000 y10=3.000000 y110=3.000000 e10=0.000000this is the second-order Runge-Kutta(middle) methodplease inpu

43、t a b and y0 ：0 1 1 please input N：10 x0=0.000000 y0=1.000000 y10=1.000000 e0=0.000000 x1=0.100000 y1=1.110000 y11=1.110000 e1=0.000000 x2=0.200000 y2=1.240000 y12=1.240000 e2=0.000000 x3=0.300000 y3=1.390000 y13=1.390000 e3=0.000000 x4=0.400000 y4=1.560000 y14=1.560000 e4=0.000000 x5=0.500000 y5=1.

44、750000 y15=1.750000 e5=0.000000 x6=0.600000 y6=1.960000 y16=1.960000 e6=0.000000 x7=0.700000 y7=2.190000 y17=2.190000 e7=0.000000 x8=0.800000 y8=2.440000 y18=2.440000 e8=0.000000 x9=0.900000 y9=2.710000 y19=2.710000 e9=0.000000 x10=1.000000 y10=3.000000 y110=3.000000 e10=0.000000（2）用三階龍格庫塔方法求解初值問題模型

45、： （4.2）結(jié)果如下：this is the third-order Heun method:please input a b and y0：0 2 1 please input N：10 x0=0.000000 y0=1.000000 y10=1.000000 e0=0.000000 x1=0.200000 y1=1.248000 y11=1.248000 e1=0.000000 x2=0.400000 y2=1.624000 y12=1.624000 e2=0.000000 x3=0.600000 y3=2.176000 y13=2.176000 e3=0.000000 x4=0.800

46、000 y4=2.952000 y14=2.952000 e4=0.000000 x5=1.000000 y5=4.000000 y15=4.000000 e5=-0.000000 x6=1.200000 y6=5.368000 y16=5.368000 e6=0.000000 x7=1.400000 y7=7.104000 y17=7.104000 e7=-0.000000 x8=1.600000 y8=9.256000 y18=9.256000 e8=-0.000000 x9=1.800000 y9=11.872000 y19=11.872000 e9=-0.000000 x10=2.00

47、0000 y10=15.000000 y110=15.000000 e10=-0.000000this is kutta methodplease input a b and y0：0 2 1 please input N：10 x0=0.000000 y0=1.000000 y10=1.000000 e0=0.000000 x1=0.200000 y1=1.248000 y11=1.248000 e1=0.000000 x2=0.400000 y2=1.624000 y12=1.624000 e2=0.000000 x3=0.600000 y3=2.176000 y13=2.176000 e

48、3=0.000000 x4=0.800000 y4=2.952000 y14=2.952000 e4=0.000000 x5=1.000000 y5=4.000000 y15=4.000000 e5=0.000000 x6=1.200000 y6=5.368000 y16=5.368001 e6=0.000000 x7=1.400000 y7=7.104000 y17=7.104000 e7=-0.000000 x8=1.600000 y8=9.256001 y18=9.256001 e8=0.000000 x9=1.800000 y9=11.872001 y19=11.872001 e9=0

49、.000000 x10=2.000000 y10=15.000001 y110=15.000000 e10=-0.0000014.3.3 結(jié)論如果常微分方程的解析解是任意的次多項式，總存在階龍格庫塔方法，用它求解該微分方程的所得的數(shù)值解與準(zhǔn)確解相符。證明：龍格庫塔方法的推導(dǎo)是基于Taylor級數(shù)方法的，階龍格庫塔方法的局部截斷誤差為：，其中：，為常數(shù)，如果常微分方程的解析解是任意的次多項式，則有，因此局部截斷誤差等于零。從而上述結(jié)論成立。5 步長算法模型及其分析下面我們分析步長對不同數(shù)值方法的影響。用歐拉方法解初值問題建立初值模型： （4.3）當(dāng)h=0.2時this is Euler meth

50、odplease input a b and y0：0 1 1 please input N： 5x0=0.000000 y0=1.000000 y10=1.000000 e0=0.000000 x1=0.200000 y1=1.200000 y11=1.183216 e1=-0.016784x2=0.400000 y2=1.373333 y12=1.341641 e2=-0.031693x3=0.600000 y3=1.531495 y13=1.483240 e3=-0.048256x4=0.800000 y4=1.681085 y14=1.612452 e4=-0.068633x5=1.0

51、00000 y5=1.826948 y15=1.732051 e5=-0.094898當(dāng)h=0.1時this is Euler method :please input a b and y0：0 1 1 please input N： 10 x0=0.000000 y0=1.000000 y10=1.000000 e0=0.000000 x1=0.100000 y1=1.100000 y11=1.095445 e1=-0.004555x2=0.200000 y2=1.191818 y12=1.183216 e2=-0.008602x3=0.300000 y3=1.277438 y13=1.26

52、4911 e3=-0.012527x4=0.400000 y4=1.358213 y14=1.341641 e4=-0.016572x5=0.500000 y5=1.435133 y15=1.414214 e5=-0.020920 x6=0.600000 y6=1.508966 y16=1.483240 e6=-0.025727x7=0.700000 y7=1.580338 y17=1.549193 e7=-0.031145x8=0.800000 y8=1.649784 y18=1.612452 e8=-0.037332x9=0.900000 y9=1.717780 y19=1.673320

53、e9=-0.044460 x10=1.000000 y10=1.784771 y110=1.732051 e10=-0.052720用改進的歐拉方法解初值問題模型（4.3）。當(dāng)h=0.2時this is improve Euler method :please input a b and y0：0 1 1 please input N：5x0=0.000000 y0=1.000000 y10=1.000000 e0=0.000000 x1=0.200000 y1=1.186667 y11=1.183216 e1=-0.003451x2=0.400000 y2=1.348312 y12=1.34

54、1641 e2=-0.006672x3=0.600000 y3=1.493704 y13=1.483240 e3=-0.010464x4=0.800000 y4=1.627861 y14=1.612452 e4=-0.015410 x5=1.000000 y5=1.754205 y15=1.732051 e5=-0.022154當(dāng)h=0.1時this is improve Euler method :please input a b and y0：0 1 1 please input N： 10 x0=0.000000 y0=1.000000 y10=1.000000 e0=0.000000

55、x1=0.100000 y1=1.095909 y11=1.095445 e1=-0.000464x2=0.200000 y2=1.184097 y12=1.183216 e2=-0.000881x3=0.300000 y3=1.266201 y13=1.264911 e3=-0.001290 x4=0.400000 y4=1.343360 y14=1.341641 e4=-0.001719x5=0.500000 y5=1.416402 y15=1.414214 e5=-0.002188x6=0.600000 y6=1.485955 y16=1.483240 e6=-0.002715x7=0.

56、700000 y7=1.552514 y17=1.549193 e7=-0.003320 x8=0.800000 y8=1.616474 y18=1.612452 e8=-0.004023x9=0.900000 y9=1.678166 y19=1.673320 e9=-0.004846x10=1.000000 y10=1.737867 y110=1.732051 e10=-0.005816用二級二階龍格庫塔方法（二階Heun方法）解初值問題模型（4.3）。當(dāng)h=0.2時this is the second-order Runge-Kutta method :please input a b a

57、nd y0：0 1 1 please input N：5x0=0.000000 y0=1.000000 y10=1.000000 e0=0.000000 x1=0.200000 y1=1.184706 y11=1.183216 e1=-0.001490 x2=0.400000 y2=1.344644 y12=1.341641 e2=-0.003003x3=0.600000 y3=1.488062 y13=1.483240 e3=-0.004822x4=0.800000 y4=1.619653 y14=1.612452 e4=-0.007201x5=1.000000 y5=1.742496 y1

58、5=1.732051 e5=-0.010445當(dāng)h=0.1時this is the second-order Runge-Kutta method :please input a b and y0：0 1 1 please input N：10 x0=0.000000 y0=1.000000 y10=1.000000 e0=0.000000 x1=0.100000 y1=1.095625 y11=1.095445 e1=-0.000180 x2=0.200000 y2=1.183572 y12=1.183216 e2=-0.000356x3=0.300000 y3=1.265449 y13=1

59、.264911 e3=-0.000538x4=0.400000 y4=1.342374 y14=1.341641 e4=-0.000733x5=0.500000 y5=1.415162 y15=1.414214 e5=-0.000948x6=0.600000 y6=1.484431 y16=1.483240 e6=-0.001191x7=0.700000 y7=1.550663 y17=1.549193 e7=-0.001470 x8=0.800000 y8=1.614246 y18=1.612452 e8=-0.001794x9=0.900000 y9=1.675493 y19=1.6733

60、20 e9=-0.002173x10=1.000000 y10=1.734671 y110=1.732051 e10=-0.002620用三級三階龍格庫塔方法（Kutta方法）解初值問題模型（4.3）。當(dāng)h=0.2時this is kutta method :please input a b and y0 ：0 1 1 please input N：5x0=0.000000 y0=1.000000 y10=1.000000 e0=0.000000 x1=0.200000 y1=1.183244 y11=1.183216 e1=-0.000028x2=0.400000 y2=1.341729 y