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

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

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

3、Student 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 differential equation solution question. Currently, the development of computer has provided the extremely powerful t

4、ool for the ordinary differential equation application and the fundamental research, the computer solving differential equation mainly uses value method. The Euler method and the RungeKutta method are the most typical commonly value method to solve the differential equation. This article dissects th

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

6、 and the difference of the coefficient different same step Rungekutta method. Through establishing C program on the computer 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 process Euler method RungeKutta method truncation error convergence目 錄微分方程常用的兩種數(shù)值解法：

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

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

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

11、們希望它產(chǎn)生的數(shù)值解收斂于初值問題的準確解，“收斂性”的一般定義為：對于所有滿足引理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)定性。當時，若舍入誤差引起的后果是有

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

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

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

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

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

17、utta-Gill)方法 (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)式容易看出，當與無關(guān)時，公式(3.3.5)實際上就是用標準的Simpson積分公式計算積分得到的，因此，可以將經(jīng)典龍格庫塔方法視為是用Simpson積分公式的推廣形式計算積分而得到

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

19、分函數(shù)，不妨設(shè)變量t為時間（以年為單位）。設(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）確定其微分方程的標準形式，這就是例4.1.1的微分方程模型。4.1.3 解決問題分別以代入，x表示利息贏利低于工資支出的數(shù)額，y表示利息贏利與工資支出平衡的數(shù)額，z表示利息贏利高于工資支出的數(shù)額，計算結(jié)果見表1.表1nn11460.25003540.13834.92625004165.0721418.4

20、25003581.614768.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.77

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

22、化，加入了計算誤差，計算結(jié)果如下（表示近似解、表示準確解、表示誤差）：This is Euler method:please input a b and y0：0 1 1 please input N：40x0=0. y0=1. y10=1. e0=0.x1=0. y1=1. y11=1. e1=-0.x2=0. y2=1. y12=1. e2=-0.x3=0. y3=1. y13=1. e3=-0.x4=0. y4=1. y14=1. e4=-0.x5=0. y5=1. y15=1. e5=-0.x6=0. y6=1. y16=1. e6=-0.x7=0. y7=1. y17=1. e7=-

23、0.x8=0. y8=1. y18=1. e8=-0.x9=0. y9=1. y19=1. e9=-0.x10=0. y10=1. y110=1. e10=-0.x11=0. y11=1. y111=1. e11=-0.x12=0. y12=1. y112=1. e12=-0.x13=0. y13=1. y113=1. e13=-0.x14=0. y14=1. y114=1. e14=-0.x15=0. y15=1. y115=1. e15=-0.x16=0. y16=1. y116=1. e16=-0.x17=0. y17=1. y117=1. e17=-0.x18=0. y18=1. y1

24、18=1. e18=-0.x19=0. y19=1. y119=1. e19=-0.x20=0. y20=1. y120=1. e20=-0.x21=0. y21=1. y121=1. e21=-0.x22=0. y22=1. y122=1. e22=-0.x23=0. y23=1. y123=1. e23=-0.x24=0. y24=1. y124=1. e24=-0.x25=0. y25=1. y125=1. e25=-0.x26=0. y26=1. y126=1. e26=-0.x27=0. y27=1. y127=1. e27=-0.x28=0. y28=1. y128=1. e28=

25、-0.x29=0. y29=1. y129=1. e29=-0.x30=0. y30=1. y130=1. e30=-0.x31=0. y31=1. y131=1. e31=-0.x32=0. y32=1. y132=1. e32=-0.x33=0. y33=1. y133=1. e33=-0.x34=0. y34=1. y134=1. e34=-0.x35=0. y35=1. y135=1. e35=-0.x36=0. y36=1. y136=1. e36=-0.x37=0. y37=1. y137=1. e37=-0.x38=0. y38=1. y138=1. e38=-0.x39=0.

26、y39=1. y139=1. e39=-0.x40=1. y40=1. y140=1. e40=-0.this is improve Euler method please input a b and y0:：0 1 1 please input N：20x0=0. y0=1. y10=1. e0=0.x1=0. y1=1. y11=1. e1=-0.x2=0. y2=1. y12=1. e2=-0.x3=0. y3=1. y13=1. e3=-0.x4=0. y4=1. y14=1. e4=-0.x5=0. y5=1. y15=1. e5=-0.x6=0. y6=1. y16=1. e6=-

27、0.x7=0. y7=1. y17=1. e7=-0.x8=0. y8=1. y18=1. e8=-0.x9=0. y9=1. y19=1. e9=-0.x10=0. y10=1. y110=1. e10=-0.x11=0. y11=1. y111=1. e11=-0.x12=0. y12=1. y112=1. e12=-0.x13=0. y13=1. y113=1. e13=-0.x14=0. y14=1. y114=1. e14=-0.x15=0. y15=1. y115=1. e15=-0.x16=0. y16=1. y116=1. e16=-0.x17=0. y17=1. y117=1

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

29、寄生蟲的增長函數(shù)如下列常微分方程組所示： 問題：預測3年后這一對寄生蟲的數(shù)量。4.2.2 模型建立分析：這是一個典型的常微分方程例題，要求預測出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預估校正公式 取h=1，計算結(jié)

30、果如表2.表2t/年u(t)v(t)11.61.221.024571.2683430.1.336640.1.41077我把龍格庫塔的C程序進行編輯后得到結(jié)果:（ei表示誤差）please input a b and y0:：0 1 1 please input N： 10x0=0. y0=1. y10=1. e0=0.x1=0. y1=1. y11=1. e1=-0.x2=0. y2=1. y12=1. e2=-0.x3=0. y3=1. y13=1. e3=-0.x4=0. y4=1. y14=1. e4=-0.x5=0. y5=1. y15=1. e5=-0.x6=0. y6=1. y16

31、=1. e6=-0.x7=0. y7=1. y17=1. e7=-0.x8=0. y8=1. y18=1. e8=-0.x9=0. y9=1. y19=1. e9=-0.x10=1. y10=1. y110=1. e10=-0.4.2.3 結(jié)論從上面的結(jié)果可以分析，用每一種方法計算節(jié)點=0.1、0.20.9、1.0，上的值都需要計算4次的值，即它們的計算量基本相同，其結(jié)果是經(jīng)典的龍格庫塔方法的精度最好，龍格庫塔方法的推導是基于Taylor級數(shù)方法的，因而在使用高階龍格庫塔方法計算時可以很精確的推算出寄生蟲每一年的數(shù)量。4.3 分別使用二階、三階龍格庫塔法解初值問題對一些特殊的微分方程，使用歐拉

32、方法和低階的龍格庫塔方法也能達到很高的精度，例如：微分方程的解析解是一次函數(shù)，則用歐拉方法求得的數(shù)值解與準確解相符，微分方程的解析解是二次函數(shù)，則用二階龍格庫塔方法求得的數(shù)值解與準確解相符。微分方程的解析解是三次多項式，用三階龍格庫塔方法求得的數(shù)值解與準確解相符。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

33、 input N：10x0=0. y0=1. y10=1. e0=0.x1=0. y1=1. y11=1. e1=0.x2=0. y2=1. y12=1. e2=0.x3=0. y3=1. y13=1. e3=0.x4=0. y4=1. y14=1. e4=0.x5=0. y5=1. y15=1. e5=0.x6=0. y6=1. y16=1. e6=0.x7=0. y7=2. y17=2. e7=0.x8=0. y8=2. y18=2. e8=0.x9=0. y9=2. y19=2. e9=0.x10=1. y10=3. y110=3. e10=0.this is improve Euler

34、 methodplease input a b and y0：0 1 1 please input N：10x0=0. y0=1. y10=1. e0=0.x1=0. y1=1. y11=1. e1=0.x2=0. y2=1. y12=1. e2=0.x3=0. y3=1. y13=1. e3=0.x4=0. y4=1. y14=1. e4=0.x5=0. y5=1. y15=1. e5=0.x6=0. y6=1. y16=1. e6=0.x7=0. y7=2. y17=2. e7=0.x8=0. y8=2. y18=2. e8=0.x9=0. y9=2. y19=2. e9=0.x10=1.

35、 y10=3. y110=3. e10=0.this is the second-order Runge-Kutta(middle) methodplease input a b and y0 ：0 1 1 please input N：10x0=0. y0=1. y10=1. e0=0.x1=0. y1=1. y11=1. e1=0.x2=0. y2=1. y12=1. e2=0.x3=0. y3=1. y13=1. e3=0.x4=0. y4=1. y14=1. e4=0.x5=0. y5=1. y15=1. e5=0.x6=0. y6=1. y16=1. e6=0.x7=0. y7=2.

36、 y17=2. e7=0.x8=0. y8=2. y18=2. e8=0.x9=0. y9=2. y19=2. e9=0.x10=1. y10=3. y110=3. e10=0.（2）用三階龍格庫塔方法求解初值問題模型： （4.2）結(jié)果如下：this is the third-order Heun method:please input a b and y0：0 2 1 please input N：10x0=0. y0=1. y10=1. e0=0.x1=0. y1=1. y11=1. e1=0.x2=0. y2=1. y12=1. e2=0.x3=0. y3=2. y13=2. e3=0.

37、x4=0. y4=2. y14=2. e4=0.x5=1. y5=4. y15=4. e5=-0.x6=1. y6=5. y16=5. e6=0.x7=1. y7=7. y17=7. e7=-0.x8=1. y8=9. y18=9. e8=-0.x9=1. y9=11. y19=11. e9=-0.x10=2. y10=15. y110=15. e10=-0.this is kutta methodplease input a b and y0：0 2 1 please input N：10x0=0. y0=1. y10=1. e0=0.x1=0. y1=1. y11=1. e1=0.x2=0

38、. y2=1. y12=1. e2=0.x3=0. y3=2. y13=2. e3=0.x4=0. y4=2. y14=2. e4=0.x5=1. y5=4. y15=4. e5=0.x6=1. y6=5. y16=5. e6=0.x7=1. y7=7. y17=7. e7=-0.x8=1. y8=9. y18=9. e8=0.x9=1. y9=11. y19=11. e9=0.x10=2. y10=15. y110=15. e10=-0.4.3.3 結(jié)論 如果常微分方程的解析解是任意的次多項式，總存在階龍格庫塔方法，用它求解該微分方程的所得的數(shù)值解與準確解相符。證明：龍格庫塔方法的推導是基于

39、Taylor級數(shù)方法的，階龍格庫塔方法的局部截斷誤差為：，其中：，為常數(shù)，如果常微分方程的解析解是任意的次多項式，則有，因此局部截斷誤差等于零。從而上述結(jié)論成立。5 步長算法模型及其分析下面我們分析步長對不同數(shù)值方法的影響。用歐拉方法解初值問題建立初值模型： （4.3）當h=0.2時this is Euler methodplease input a b and y0：0 1 1 please input N： 5x0=0. y0=1. y10=1. e0=0.x1=0. y1=1. y11=1. e1=-0.x2=0. y2=1. y12=1. e2=-0.x3=0. y3=1. y13=1

40、. e3=-0.x4=0. y4=1. y14=1. e4=-0.x5=1. y5=1. y15=1. e5=-0.當h=0.1時this is Euler method :please input a b and y0：0 1 1 please input N： 10x0=0. y0=1. y10=1. e0=0.x1=0. y1=1. y11=1. e1=-0.x2=0. y2=1. y12=1. e2=-0.x3=0. y3=1. y13=1. e3=-0.x4=0. y4=1. y14=1. e4=-0.x5=0. y5=1. y15=1. e5=-0.x6=0. y6=1. y16=

41、1. e6=-0.x7=0. y7=1. y17=1. e7=-0.x8=0. y8=1. y18=1. e8=-0.x9=0. y9=1. y19=1. e9=-0.x10=1. y10=1. y110=1. e10=-0.用改進的歐拉方法解初值問題模型（4.3）。當h=0.2時this is improve Euler method :please input a b and y0：0 1 1 please input N：5x0=0. y0=1. y10=1. e0=0.x1=0. y1=1. y11=1. e1=-0.x2=0. y2=1. y12=1. e2=-0.x3=0. y3=

42、1. y13=1. e3=-0.x4=0. y4=1. y14=1. e4=-0.x5=1. y5=1. y15=1. e5=-0.當h=0.1時this is improve Euler method :please input a b and y0：0 1 1 please input N： 10x0=0. y0=1. y10=1. e0=0.x1=0. y1=1. y11=1. e1=-0.x2=0. y2=1. y12=1. e2=-0.x3=0. y3=1. y13=1. e3=-0.x4=0. y4=1. y14=1. e4=-0.x5=0. y5=1. y15=1. e5=-0.

43、x6=0. y6=1. y16=1. e6=-0.x7=0. y7=1. y17=1. e7=-0.x8=0. y8=1. y18=1. e8=-0.x9=0. y9=1. y19=1. e9=-0.x10=1. y10=1. y110=1. e10=-0.用二級二階龍格庫塔方法（二階Heun方法）解初值問題模型（4.3）。當h=0.2時this is the second-order Runge-Kutta method :please input a b and y0：0 1 1 please input N：5x0=0. y0=1. y10=1. e0=0.x1=0. y1=1. y11

44、=1. e1=-0.x2=0. y2=1. y12=1. e2=-0.x3=0. y3=1. y13=1. e3=-0.x4=0. y4=1. y14=1. e4=-0.x5=1. y5=1. y15=1. e5=-0.當h=0.1時this is the second-order Runge-Kutta method :please input a b and y0：0 1 1 please input N：10x0=0. y0=1. y10=1. e0=0.x1=0. y1=1. y11=1. e1=-0.x2=0. y2=1. y12=1. e2=-0.x3=0. y3=1. y13=1

45、. e3=-0.x4=0. y4=1. y14=1. e4=-0.x5=0. y5=1. y15=1. e5=-0.x6=0. y6=1. y16=1. e6=-0.x7=0. y7=1. y17=1. e7=-0.x8=0. y8=1. y18=1. e8=-0.x9=0. y9=1. y19=1. e9=-0.x10=1. y10=1. y110=1. e10=-0.用三級三階龍格庫塔方法（Kutta方法）解初值問題模型（4.3）。當h=0.2時this is kutta method :please input a b and y0 ：0 1 1 please input N：5x0=0. y0=1. y10=1. e0=0.x1=0. y1=1. y11=1. e1=-0.x2=0. y2=1. y12=1. e2=-0.x3=0. y3=1. y13=1. e3=-0.x4=0. y4=1.