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

2、它們的優(yōu)缺點(diǎn)，以及它們的收斂性，相容性，及穩(wěn)定性。討論了步長(zhǎng)的變化對(duì)數(shù)值方法的影響和系數(shù)不同的同階龍格庫(kù)塔方法的差別。通過編制C程序在計(jì)算機(jī)上實(shí)現(xiàn)這兩類方法及對(duì)一些典型算例的結(jié)果分析比較，能更深切體會(huì)它們的功能，優(yōu)缺點(diǎn)及適用場(chǎng)合，從而在實(shí)際應(yīng)用中能對(duì)不同類型和不同要求的常微分方程會(huì)選取適當(dāng)?shù)那蠼夥椒āｊP(guān)鍵詞：顯式單步法 歐拉（Euler）方法 龍格庫(kù)塔（RungeKutta）方法 截?cái)嗾`差 收斂性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、歐拉方法與龍格庫(kù)塔法前言 常微分方程的形成與發(fā)展是和力學(xué)、天文學(xué)、物理學(xué)以及其他科學(xué)技術(shù)的發(fā)展密切相關(guān)的。數(shù)學(xué)其他分支的新發(fā)展，如復(fù)變函數(shù)、群、組合拓?fù)鋵W(xué)等，都對(duì)常微分方程的發(fā)展產(chǎn)生了深刻的影響，當(dāng)前計(jì)算機(jī)的發(fā)展更是為常微分方程的應(yīng)用及理論研究提供了非常有力的工具。牛頓在研究天體力學(xué)和機(jī)械力學(xué)的時(shí)候，利用了微分方程這個(gè)工具，從理論上得到了行星運(yùn)動(dòng)規(guī)律。后來(lái)，法國(guó)天文學(xué)家勒維烈和英國(guó)天文學(xué)家亞當(dāng)斯使用微分方程各自計(jì)算出那時(shí)尚未發(fā)現(xiàn)的海王星的位置。這些都使數(shù)學(xué)家更加深信微分方程在認(rèn)識(shí)自然、改造自然方面的巨大力量，微分方程也就成了最有生命力的數(shù)學(xué)分支。然而，我們知道，只有少數(shù)十分簡(jiǎn)單的微分方程能夠用

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

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

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

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

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

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

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

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

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

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

19、分函數(shù)，不妨設(shè)變量t為時(shí)間（以年為單位）。設(shè)為第n年的凈資產(chǎn)，為第n+1年的凈資產(chǎn)，利息以每年4%的速度增長(zhǎng)，且公司每年支付職工工資為100萬(wàn)，則第n+1年的凈資產(chǎn)增長(zhǎng)數(shù)額為，由于已得第n年的凈資產(chǎn)為，則第n年的凈資產(chǎn)加上第n+1年的凈資產(chǎn)增長(zhǎng)數(shù)額就得到第n+1年的凈資產(chǎn)。歸納公式： （4-1）確定其微分方程的標(biāo)準(zhǔn)形式，這就是例4.1.1的微分方程模型。4.1.3 解決問題分別以代入，x表示利息贏利低于工資支出的數(shù)額，y表示利息贏利與工資支出平衡的數(shù)額，z表示利息贏利高于工資支出的數(shù)額，計(jì)算結(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可以看到當(dāng)利息贏利低于工資的支出，公司的凈資產(chǎn)逐年減少，以致凈資產(chǎn)為負(fù)值；當(dāng)利息贏利與工資的支出平衡時(shí)，公司的凈資產(chǎn)每年保持不變；當(dāng)利息贏利超過工資的支出，公司的凈資產(chǎn)穩(wěn)步增長(zhǎng)。再用歐拉方法求解，每計(jì)算一步，依次需要計(jì)算1次、2次和4次函數(shù)的值，為了比較在計(jì)算量相同的條件下近似解的精度，步長(zhǎng)分別取0.025、0.05、0.1。我對(duì)歐拉的C程序做了一些優(yōu)

22、化，加入了計(jì)算誤差，計(jì)算結(jié)果如下（表示近似解、表示準(zhǔn)確解、表示誤差）：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)鍵在于如何預(yù)測(cè)若干年后的凈資產(chǎn)趨勢(shì)，用普通的微分方程就無(wú)法進(jìn)行預(yù)測(cè)，且人工計(jì)算量相當(dāng)大，這里我使用歐拉方法可以計(jì)算出精度以及誤差，通過電腦運(yùn)行程序就可以預(yù)測(cè)出若干年后的凈資產(chǎn)情況，歐拉方法的使用使得解題更加方便且精確。4.2 龍格庫(kù)塔解題及其數(shù)學(xué)模型4.2.1 問題提出 例4.2.1 兩種果樹寄生蟲，其數(shù)量分別是其中一種寄生蟲以吃另一種寄生蟲為生，兩種

29、寄生蟲的增長(zhǎng)函數(shù)如下列常微分方程組所示： 問題：預(yù)測(cè)3年后這一對(duì)寄生蟲的數(shù)量。4.2.2 模型建立分析：這是一個(gè)典型的常微分方程例題，要求預(yù)測(cè)出3年后這對(duì)寄生蟲的數(shù)量。確定變量：假定時(shí)間是寄生蟲數(shù)量的積分函數(shù)，不妨設(shè)變量t為時(shí)間（以年為單位）。由題知有兩種寄生蟲u和v，u寄生蟲年增長(zhǎng)函數(shù)為，v寄生蟲年增長(zhǎng)函數(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，計(jì)算結(jié)

30、果如表2.表2t/年u(t)v(t)11.61.221.024571.2683430.1.336640.1.41077我把龍格庫(kù)塔的C程序進(jìn)行編輯后得到結(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é)果可以分析，用每一種方法計(jì)算節(jié)點(diǎn)=0.1、0.20.9、1.0，上的值都需要計(jì)算4次的值，即它們的計(jì)算量基本相同，其結(jié)果是經(jīng)典的龍格庫(kù)塔方法的精度最好，龍格庫(kù)塔方法的推導(dǎo)是基于Taylor級(jí)數(shù)方法的，因而在使用高階龍格庫(kù)塔方法計(jì)算時(shí)可以很精確的推算出寄生蟲每一年的數(shù)量。4.3 分別使用二階、三階龍格庫(kù)塔法解初值問題對(duì)一些特殊的微分方程，使用歐拉

32、方法和低階的龍格庫(kù)塔方法也能達(dá)到很高的精度，例如：微分方程的解析解是一次函數(shù)，則用歐拉方法求得的數(shù)值解與準(zhǔn)確解相符，微分方程的解析解是二次函數(shù)，則用二階龍格庫(kù)塔方法求得的數(shù)值解與準(zhǔn)確解相符。微分方程的解析解是三次多項(xiàng)式，用三階龍格庫(kù)塔方法求得的數(shù)值解與準(zhǔn)確解相符。4.3.1 建立模型建立初值問題模型： (4.1)4.3.2 用不同的龍格庫(kù)塔法解（4.1）初值問題模型 （1）用二階龍格庫(kù)塔方法求解初值問題（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）用三階龍格庫(kù)塔方法求解初值問題模型： （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é)論 如果常微分方程的解析解是任意的次多項(xiàng)式，總存在階龍格庫(kù)塔方法，用它求解該微分方程的所得的數(shù)值解與準(zhǔn)確解相符。證明：龍格庫(kù)塔方法的推導(dǎo)是基于

39、Taylor級(jí)數(shù)方法的，階龍格庫(kù)塔方法的局部截?cái)嗾`差為：，其中：，為常數(shù)，如果常微分方程的解析解是任意的次多項(xiàng)式，則有，因此局部截?cái)嗾`差等于零。從而上述結(jié)論成立。5 步長(zhǎng)算法模型及其分析下面我們分析步長(zhǎng)對(duì)不同數(shù)值方法的影響。用歐拉方法解初值問題建立初值模型： （4.3）當(dāng)h=0.2時(shí)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.當(dāng)h=0.1時(shí)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.用改進(jìn)的歐拉方法解初值問題模型（4.3）。當(dāng)h=0.2時(shí)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.當(dāng)h=0.1時(shí)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.用二級(jí)二階龍格庫(kù)塔方法（二階Heun方法）解初值問題模型（4.3）。當(dāng)h=0.2時(shí)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.當(dāng)h=0.1時(shí)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.用三級(jí)三階龍格庫(kù)塔方法（Kutta方法）解初值問題模型（4.3）。當(dāng)h=0.2時(shí)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.