一階微分方程的matlab數(shù)值解法

1、*大學(xué)畢業(yè)設(shè)計(jì)-PAGE - 2 -. z. - - .可修編 .畢業(yè)論文題 目一階微分方程的matlab數(shù)值解法學(xué) 院數(shù)學(xué)科學(xué)學(xué)院專 業(yè)數(shù)學(xué)與應(yīng)用數(shù)學(xué)班 級(jí)數(shù)學(xué)1303班 學(xué) 生鄒健峰學(xué) 號(hào)7指導(dǎo)教師邢順來(lái)講師二一七年五月七日-PAGE . z.摘 要微分方程的數(shù)值解法是近現(xiàn)代數(shù)學(xué)家和科學(xué)家們研究的熱點(diǎn)，微分方程的MATLAB數(shù)值解法可以幫助我們解決現(xiàn)實(shí)生活中的許多數(shù)學(xué)問(wèn)題，提高計(jì)算機(jī)幫助我們解決問(wèn)題的效率。本文主要討論研究一階微分方程的MATLAB數(shù)值解法中的三種Euler法和三種Runge-Kutta法，介紹以上六種方法的主要容，簡(jiǎn)單介紹了MATLAB的相關(guān)容和微分方程的開(kāi)展簡(jiǎn)史。通過(guò)具

2、體的微分方程來(lái)研究以上算法的編程實(shí)現(xiàn)，通過(guò)MATBLAB求解具體的一階微分方程來(lái)探究以上方法的優(yōu)缺點(diǎn)，通過(guò)圖表來(lái)分析得出結(jié)論：改進(jìn)Euler法和四階Runge-Kutta法的階的精度較高，具有較好的算法穩(wěn)定性。關(guān)鍵詞：一階微分方程；數(shù)值解法；MATLAB；Euler法；Runge-Kutta法ABSTRACTThe numerical solution of differential equations is a hot spot for modern mathematicians and scientists,the MATLAB numerical solution of the diff

3、erential equation can help us solve many mathematical problems in real life,improve the efficiency of the puter to help us solve the problem.This paper mainly discusses three kinds of Euler methods and three Runge-Kutta methods in MATLAB numerical solution of first order differential equation,introd

4、uce the main contents of the above si* methods, a brief introduction to the development of MATLAB related content and differential equations of a brief history.The programming of the above algorithm is studied by the concrete differential equation,he advantages and disadvantages of the above methods

5、 are e*plored through MATBLAB solving the specific first-order differential equation,through the chart to analyze the conclusions:The improved Euler method and the fourth order Runge-Kutta method have higher accuracy，has good algorithm stability.Keywords：First Order Differential Equation; Numerical

6、Solution; MATLAB; Euler Method; Runge-Kutta Method目錄摘要 PAGEREF _Toc29420 IABSTRACT PAGEREF _Toc434 II1 前言 PAGEREF _Toc19134 11.1 引言 PAGEREF _Toc14582 11.2常微分方程的開(kāi)展簡(jiǎn)史 PAGEREF _Toc21991 21.3 國(guó)外研究現(xiàn)狀 PAGEREF _Toc15673 21.4 研究主要容及研究方案 PAGEREF _Toc11775 31.4.1 研究的主要容 PAGEREF _Toc1002 3 研究方案 PAGEREF _Toc297

7、35 31.5 研究難點(diǎn)32 預(yù)備知識(shí)42.1顯式Euler法42.2隱式Euler法52.3改進(jìn)Euler法52.4二階Runge-Kutta法62.5三階Runge-Kutta法72.6四階Runge-Kutta法72.7單步法的收斂性和穩(wěn)定性93 微分方程數(shù)值解法的編程實(shí)現(xiàn)113.1 常微分方程的符號(hào)解法113.2顯式Euler法的編程實(shí)現(xiàn)123.3隱式Euler法的編程實(shí)現(xiàn)143.4改進(jìn)Euler法的編程實(shí)現(xiàn)173.5二階Runge-Kutta法的編程實(shí)現(xiàn)203.6三階Runge-Kutta法的編程實(shí)現(xiàn)223.7四階Runge-Kutta法的編程實(shí)現(xiàn)24結(jié)論29參考文獻(xiàn)30致31附錄3

8、2-. z.1前言1.1引言James Bernoulli在1676年寫給Newton的信中首次談到了微分方程。大約到了十八世紀(jì)中期，微分方程才成為一門獨(dú)立的學(xué)科，微分方程是研究和提醒自然規(guī)律不可或缺的重要工具。眾所周知，在1846年數(shù)學(xué)家和天文學(xué)家共同發(fā)現(xiàn)的海王星就是微分方程的功績(jī)。還有在1991年，科學(xué)考察隊(duì)在Alps山脈發(fā)現(xiàn)了一個(gè)冰封保存完好的人類，然后通過(guò)儀器測(cè)得碳元素衰變的速率，建立微分方程數(shù)學(xué)模型，得出了這個(gè)人是生活在距今5000多年前的時(shí)代。Newton、Leibniz、Bernoulli family、Lagrange、Laplace等著名數(shù)學(xué)家們都對(duì)微分方程的開(kāi)展做出了巨大的

9、奉獻(xiàn)。微分方程的定義：如果知道自變量、未知函數(shù)以及函數(shù)的導(dǎo)數(shù)或微分組成的關(guān)系式我們就稱為微分方程。常微分方程的定義：通過(guò)求解微分方程得出未知函數(shù)，并且自變量只有一個(gè)的微分方程我們就叫做常微分方程，自變量個(gè)數(shù)為兩個(gè)或者兩個(gè)以上的微分方程稱為偏微分方程1。我們把含有n個(gè)互不相等的任意常數(shù)的解稱為n階常微分方程的通解。微分方程的特解是指滿足微分方程初值條件的解1。我們?cè)谏钪袝?huì)經(jīng)常用到一些常微分方程模型，例如英國(guó)著名人口統(tǒng)計(jì)學(xué)家Malthus提出的著名的Malthus人口模型；荷蘭著名生物學(xué)家Verhulst提出的logistic人口增長(zhǎng)模型；傳染病模型SI模型、SIS模型無(wú)免疫性、SIR強(qiáng)免疫性模

10、型；兩生物種群生態(tài)模型；氣象學(xué)家Lorenz提出的Lorenz方程用來(lái)預(yù)測(cè)天氣變化；化學(xué)動(dòng)力學(xué)模型2?？梢?jiàn)微分方程與我們的生活有著密不可分的關(guān)系。MATLAB是分別取了Matri*和Laboratory兩個(gè)單詞的前三個(gè)字母組成的,是矩陣實(shí)驗(yàn)室的名稱縮寫。MATLAB是來(lái)自美國(guó)大學(xué)的計(jì)算機(jī)教師Cleve Moler創(chuàng)造的，并且由Cleve Moler和JohnLillte把其推向了市場(chǎng)，并得到了廣闊用戶的喜愛(ài)。MATLAB是數(shù)學(xué)科學(xué)和工程科學(xué)的高級(jí)計(jì)算機(jī)語(yǔ)言，用戶可以用數(shù)學(xué)語(yǔ)言來(lái)編寫程序，編程邏輯非常類似人們?nèi)粘鴮憯?shù)學(xué)公式的邏輯，具有簡(jiǎn)單性、廣泛性、開(kāi)放性、效率高、價(jià)格低等特點(diǎn)。并且MATLA

11、B工具箱儲(chǔ)存的MATLAB函數(shù)庫(kù)里面有大量的函數(shù)程序，可以解決許多的數(shù)學(xué)和工程科學(xué)問(wèn)題。從2010版開(kāi)場(chǎng)MATLAB的功能越來(lái)越多，編譯器可以把M文件轉(zhuǎn)換為C語(yǔ)言文件或者其他幾種常用的計(jì)算機(jī)文件，并給用戶提供源程序代碼進(jìn)展自我開(kāi)發(fā)，結(jié)合Math Works公司提供的數(shù)據(jù)庫(kù)，用戶可以利用MATLAB來(lái)開(kāi)發(fā)出功能強(qiáng)大的程序3,4,5。微分方程的MATLAB數(shù)值解法可以幫助我們解決許多現(xiàn)實(shí)生活中的數(shù)學(xué)問(wèn)題，提高計(jì)算機(jī)幫助我們解決問(wèn)題的效率，如今微分方程的MATLAB數(shù)值解法已經(jīng)開(kāi)展得比較成熟，有許多關(guān)于數(shù)學(xué)領(lǐng)域的MATLAB書籍可供我們選擇閱讀，很好的解決了我們的學(xué)習(xí)需求和使用問(wèn)題。MATLAB應(yīng)用

12、領(lǐng)域也非常廣泛，比方在數(shù)學(xué)建模和數(shù)學(xué)實(shí)驗(yàn)領(lǐng)域，我們可以通過(guò)直接調(diào)用MATLAB里面的微分方程函數(shù)文件，在操作界面輸入初始數(shù)據(jù)，還把抽象復(fù)雜的數(shù)學(xué)問(wèn)題變成簡(jiǎn)單形象的二維或者三維圖形，來(lái)幫助我們分析問(wèn)題，找出解題方法，從而解決我們生活中負(fù)責(zé)的數(shù)學(xué)問(wèn)題，表達(dá)了數(shù)學(xué)源自生活又聯(lián)系生活的主要特征。除此之外，在物理學(xué)、工程科學(xué)等科學(xué)領(lǐng)域中MATLAB也具有重要地位6,7。1.2常微分方程的開(kāi)展簡(jiǎn)史常微分方程的開(kāi)展可以簡(jiǎn)單的分為以下幾個(gè)階段：“求通解階段；“求定解階段；“求所有解階段；“求特殊解階段。在常微分方程的“求通解階段，人們希望能夠用初等函數(shù)來(lái)表示解，Newton就希望能夠用別離變量法來(lái)求解一階微分

13、方程問(wèn)題，但Euler卻嘗試用積分因子來(lái)求解一階微分方程問(wèn)題。而B(niǎo)ernoulli在研究微分方程的過(guò)程中提出了著名的Bernoulli方程，Riccati也在研究過(guò)程中提出了Riccati方程，他們的發(fā)現(xiàn)對(duì)微分方程未來(lái)的開(kāi)展打下了根底，起到了積極的作用。直到Liouville證明了Riccati方程沒(méi)有一般的初等解，還有Cauchy提出了微分方程的初值問(wèn)題，常微分方程進(jìn)入了“求定解階段。直到19世紀(jì)末，由于天體力學(xué)的研究需求，常微分方程從“求定解階段進(jìn)入了“求所有解的階段。首先是Poincar提出了定性理論，接著是Hilbert提出了關(guān)于極限環(huán)個(gè)數(shù)的問(wèn)題，促進(jìn)了微分方程定性理論的開(kāi)展，緊接著L

14、yapunov提出了運(yùn)動(dòng)穩(wěn)定性理論，廣泛應(yīng)用在物理和工程學(xué)中。而在動(dòng)力系統(tǒng)方面，Birkhoff、Arnold、Smale都做出了奉獻(xiàn)。到了20世紀(jì)70年代計(jì)算機(jī)技術(shù)的開(kāi)展帶動(dòng)了常微分方程的開(kāi)展，從“求所有解階段進(jìn)入到了“求特定解階段，發(fā)現(xiàn)了一些特定解和方程，比方混沌解、孤立子、奇異吸引子、分形。關(guān)于常微分方程的研究還在繼續(xù)，如今研究的領(lǐng)域更加廣泛，更多的微分方程模型被提出和應(yīng)用起來(lái)了1,8,9,10。 1.3國(guó)外研究現(xiàn)狀關(guān)于微分方程數(shù)值解法的研究，在沒(méi)有MATLAB以及類似的計(jì)算機(jī)軟件問(wèn)世之前，數(shù)學(xué)家和科學(xué)家們只能研究和解決一些低階的微分方程問(wèn)題，并且有些方程的解題過(guò)程過(guò)于繁瑣和復(fù)雜。但是M

15、ATLAB的出現(xiàn)促進(jìn)了數(shù)學(xué)科學(xué)研究的開(kāi)展，對(duì)數(shù)學(xué)和相關(guān)領(lǐng)域的開(kāi)展做出了巨大奉獻(xiàn)。微分方程的MATLAB數(shù)值解法也迅速開(kāi)展起來(lái)，國(guó)外現(xiàn)在主要研究通過(guò)MATLAB等計(jì)算機(jī)軟件來(lái)研究：一些高階微分方程的解的有關(guān)問(wèn)題；一些復(fù)雜的微分方程是否存在定解問(wèn)題；利用MATLAB來(lái)計(jì)算一些著名的數(shù)學(xué)模型問(wèn)題，比方Volterra被捕食 捕食模型、競(jìng)爭(zhēng)模型、共生模型、商品供需模型等等；研究混沌解、孤立子。在以上幾方面的研究上都取得了不小的進(jìn)步，使常微分方程的體系更加的完善了。在國(guó)，2009年中國(guó)數(shù)學(xué)會(huì)成功舉辦了第四屆全國(guó)青年常微分方程學(xué)術(shù)會(huì)議，許多國(guó)知名學(xué)者和學(xué)術(shù)研究者都到場(chǎng)參加學(xué)術(shù)交流，這一項(xiàng)活動(dòng)促進(jìn)了我國(guó)常微

16、分方程領(lǐng)域的開(kāi)展，為青年學(xué)者提供了珍貴經(jīng)歷，這次會(huì)議也受到了國(guó)外數(shù)學(xué)愛(ài)好者的高度關(guān)注。在書籍方面，我國(guó)也發(fā)行了許多關(guān)于微分方程數(shù)值分析和微分方程數(shù)值解法方面的書籍，都有詳細(xì)介紹微分方程關(guān)于MATLAB的數(shù)值解法，比方林成森編著的第二版?數(shù)值計(jì)算方法?里面就有詳細(xì)介紹常微分方程初值問(wèn)題的數(shù)值解法，包括有：離散變量法和離散誤差法；單步法；單步法的相容性、收斂性和穩(wěn)定性；多步法；差分方程的簡(jiǎn)介；線性步法的相容性、收斂性和數(shù)值穩(wěn)定性；常微分方程組和高階微分方程的數(shù)值解法等幾個(gè)方面，成為了21世紀(jì)高等院校的教材。周品、何正風(fēng)等編著的?MATLAB數(shù)值分析?作為了理工科各專業(yè)本科生、研究生以及應(yīng)用MATL

17、AB相關(guān)科技人員學(xué)習(xí)MATLAB數(shù)值分析、建模、仿真的教材或者參考書?，F(xiàn)如今，由于信息時(shí)代的飛速開(kāi)展，數(shù)值計(jì)算方面的計(jì)算機(jī)軟件越來(lái)越多，功能也越來(lái)越全面，所以微分方程的數(shù)值解法受到了更多人的重視，開(kāi)展前景變得更好了，這對(duì)于學(xué)習(xí)微分方程的數(shù)值解法的人來(lái)說(shuō)是非常好的消息。1.4研究的主要容與方法1.4.1 研究的主要容本文針對(duì)一階常微分方程的數(shù)值解法中的MATLAB數(shù)值解法進(jìn)展算法介紹和編程實(shí)現(xiàn)，利用MATLAB對(duì)一階微分方程中的Euler法、隱式Euler法、改進(jìn)Euler法以及Runge-Kutta法進(jìn)展精度分析并比較，提醒以上幾種方法的優(yōu)缺點(diǎn)，編程實(shí)現(xiàn)以上方法的數(shù)值解法，并通過(guò)具體事例來(lái)驗(yàn)證

18、程序的完整性、穩(wěn)定性、以及缺乏。1.4.2研究方法一階微分方程的MATLAB數(shù)值解法的編程實(shí)現(xiàn)需要通過(guò)Euler法、隱式Euler法、改進(jìn)Euler法以及二到四階Runge-Kutta法等方法的算法邏輯，針對(duì)不同一階微分方程來(lái)編寫MATLAB程序，然后再反復(fù)調(diào)試程序，找出其中的缺乏并改進(jìn)，然后保存程序并做記錄，以便下次使用時(shí)直接調(diào)用。1.5研究難點(diǎn)1微分方程的模型多樣，研究方向多，不同的模型有不同的數(shù)值分析和數(shù)值解法；2微分方程的MATLAB數(shù)值解法需要符合實(shí)際問(wèn)題，限制條件比較多；3微分方程的類型多樣，MATLAB的函數(shù)庫(kù)里只有局部微分方程的調(diào)用函數(shù)，不能滿足研究的需求，需要自定義調(diào)用函數(shù)；

19、4微分方程的MATLAB數(shù)值解法需要反復(fù)調(diào)試，得出滿意結(jié)果花費(fèi)時(shí)間較長(zhǎng)。2預(yù)備知識(shí)2.1顯式Euler法顯示Euler法也被人們簡(jiǎn)稱為Euler法，是解決微分方程初值問(wèn)題最簡(jiǎn)單的數(shù)值方法。初值問(wèn)題的解代表通過(guò)點(diǎn)的一條稱為微分方程的積分曲線。積分曲線上每一點(diǎn)的切線斜率等于函數(shù)在這點(diǎn)的值。Euler方法的求解過(guò)程是從初始點(diǎn)出發(fā),作積分曲線在點(diǎn)的切線其斜率為，與直線相交與點(diǎn)，得到作為的近似值，過(guò)點(diǎn)，以為斜率的切線方程為當(dāng)時(shí)，得這樣就獲得了點(diǎn)的坐標(biāo)3,11。不斷重復(fù)上述過(guò)程，就可以得到一系列的點(diǎn)。對(duì)已求得的點(diǎn)以為斜率做直線，當(dāng)時(shí)，便可得到方程取，這樣從逐一算出所對(duì)應(yīng)的數(shù)值解，同樣也獲得了一條曲線作為原

20、始曲線的近似曲線，因此又稱Euler法為折線法。在計(jì)算時(shí)，一般取常數(shù)，則顯式Euler公式為.2.12.2隱式Euler法假設(shè)我們用近似等式來(lái)替代顯式Euler公式右端的積分項(xiàng)，就得到再用替代，替代，我們就可以得到這就是隱式Euler公式。后退的Euler法也被稱為隱式Euler法。而我們經(jīng)常用到的梯形公式。對(duì)方程的兩端在區(qū)間上積分，得要通過(guò)這個(gè)積分關(guān)系式得的近似值，只要近似的算出其中的積分項(xiàng)。假設(shè)使用梯形方法計(jì)算積分項(xiàng)再代入中，得將式子中的分別用代替，有.2.2這種差分格式就稱為梯形格式，梯形格式上是顯示Euler公式與隱式Euler公式的算術(shù)平均。梯形方法是隱式單步法，可以用迭代法求解，與

21、后退的Euler法一樣，仍用Euler法進(jìn)展初值迭代，所以梯形公式的迭代公式為2.32.3改進(jìn)Euler法顯示Euler公式計(jì)算工作量小，但精度不高，假設(shè)對(duì)計(jì)算結(jié)果的精度要求高，就必須使用就是精度更高的計(jì)算方法。微分方程初值問(wèn)題的梯形公式相比顯示Euler公式，精度提高了，但是用迭代的算法來(lái)求解計(jì)算工作量大。所以可以改進(jìn)思路：綜合表示Euler公式和梯形公式便可以得到改進(jìn)的Euler公式。先用顯示Euler公式求出一個(gè)初步的近似值，記為，稱之為預(yù)測(cè)值，它的精度不高，再用梯形公式對(duì)它進(jìn)展一次矯正，也就是迭代一次，求出，稱為矯正值12,13。這種預(yù)測(cè)矯正的方法我們就叫做改進(jìn)的Euler法.2.4它

22、也可以表示為嵌套模式.或者表示為以下的平均化形式2.52.4二階Runge-Kutta法二階Runge-Kutta的一般公式為2.6局部截?cái)嗾`差為.其中加權(quán)系數(shù)和選點(diǎn)調(diào)節(jié)系數(shù)都是特定系數(shù)。本質(zhì)是用差商代替導(dǎo)數(shù)的情況下，把在區(qū)間上的假設(shè)干點(diǎn)的值加權(quán)平均，再作為斜率，作過(guò)點(diǎn)的直線交直線，求得下一點(diǎn)。適中選擇是提高精度的方法14,15,16。我們常用的二階Runge-Kutta公式有中點(diǎn)方法2.7和Heun方法2.82.5 三階Runge-Kutta法三階Runge-Kutta一般公式為2.9局部截?cái)嗾`差.常用的三階Runge-Kutta公式有2.10和2.112.6 四階Runge-Kutta法繼

23、續(xù)以上的演算過(guò)程，可以推導(dǎo)出各種四階Runge-Kutta公式，下面介紹一個(gè)被被稱為經(jīng)典Runge-Kutta公式，其公式的格式為2.12其截?cái)嗾`差為.四階Runge-Kutta法要求所求的解具有比較好的光滑性，反之，使用四階Runge-Kutta法求得的解精度比改進(jìn)的Euler法求得的解精度低。2.7單步法的收斂性和穩(wěn)定性定義1假設(shè)微分方程的右端函數(shù)在帶形區(qū)域中連續(xù)，并且關(guān)于滿足Lipschitz條件.假設(shè)對(duì)所有的，則稱該單步法是收斂的.定義2如果存在正常數(shù)以及，使得對(duì)任意的初始出發(fā)值，*單步法的相應(yīng)準(zhǔn)確解以及，對(duì)所有的，恒有則稱該單步法是穩(wěn)定的.定義3對(duì)給定的微分方程和給定的步長(zhǎng)h，如果由

24、單步法顯式或隱式計(jì)算時(shí)有大小為的誤差，即，而引起其后值的變化小于，則說(shuō)該單步法是絕對(duì)穩(wěn)定的.一般只限于典型的微分方程考慮數(shù)值方法的絕對(duì)穩(wěn)定性我們僅限于為實(shí)數(shù)的情形.假設(shè)對(duì)所有，單步法都絕對(duì)穩(wěn)定，則稱為絕對(duì)穩(wěn)定區(qū)間3.3微分方程數(shù)值解法的編程實(shí)現(xiàn)3.1 常微分方程的符號(hào)解法在我們開(kāi)場(chǎng)著手編寫程序前，我們先一起了解一下什么是算法，算法是計(jì)算機(jī)編程的基石。算法的定義是：算法是指解題方案的準(zhǔn)確而完整的描述，是一系列解決問(wèn)題的清晰指令，算法是一種用系統(tǒng)的方法描述如何解決問(wèn)題的策略機(jī)制3,14,15。比方Euler算法向我們介紹了如何運(yùn)用Euler法來(lái)求解問(wèn)題，是一種解決問(wèn)題的思路，只有先掌握了算法邏輯，

25、我們才可以根據(jù)具體問(wèn)題需求來(lái)編寫Euler法的MATLAB程序，或者編寫其他編程軟件的程序。而且根據(jù)實(shí)際問(wèn)題的不同，我們可以靈活的使用不同的算法。我們接著來(lái)了解一下什么是常微分方程的符號(hào)解法。當(dāng)一個(gè)微分方程能夠求得它的解析解通解時(shí)，我們就可以通過(guò)MATLAB工具箱里面的函數(shù)求出該微分方程的準(zhǔn)確解。但在實(shí)際生活中我們更需要求得微分方程的近似解或者特解，我們也可以通過(guò)MATLB求出它的數(shù)值解。一階常微分方程First-order Ordinary Differential Equation, ODE寫為其解為該解滿足，并在初始條件下，可求得唯一解。在MATLAB中，常微分方程的符號(hào)解可使用dsol

26、ve函數(shù)求得，調(diào)用格式為dsolve(eqns)； %求出變量為t的微分方程的解析解dsolve(eqns,conds)； %求出變量為t的微分方程初值問(wèn)題的解析解dsolve(eqns,Value)； %求出變量為*的微分方程的解析解dsolve(eqns,conds,Name,Value)； %求出變量為*的微分方程初值問(wèn)題的解析解在MATLAB程序代碼中，D表示對(duì)獨(dú)立變量的微分，即，所以用戶不能再自定義包含D的變量。微分方程在輸入時(shí)，要輸入成“Dy，要輸入成“D2y假設(shè)初始條件的數(shù)目比被微量數(shù)目少，則結(jié)果中將包括不定常數(shù)等。系統(tǒng)默認(rèn)“t為自變量，如果所求的方程中自變量是“*，必須要在語(yǔ)句

27、的末尾定義該方程中的自變量為“*，否則系統(tǒng)會(huì)把“t認(rèn)為自變量，從而造成錯(cuò)誤，所以要引起注意。我們通過(guò)兩個(gè)例子來(lái)了解一下dsolve函數(shù)的使用方法：1、計(jì)算微分方程的解在MATLAB命令窗口輸入clear allz=dsolve(Dy=-*,D*=y)運(yùn)行程序輸出結(jié)果z = *: 1*1 sym y: 1*1 sym z.* %查詢*值ans =C2*cos(t) - C1*sin(t) z.y %查詢y值ans = C1*cos(t) - C2*sin(t)2、計(jì)算微分方程的通解在命令窗口輸入以下命令clear ally=dsolve(Dy=*2*y,*)輸出結(jié)果y = C2*e*p(*2)所

28、以，通解為3.2Euler法的編程實(shí)現(xiàn)Euler法算法見(jiàn)附錄13、用Euler法求解一階常微分方程的數(shù)值解編寫Euler法的MATLAB程序代碼如下：function T= Euler(f,*0,y0,*n,N)% Euler.m函數(shù)為用Euler法求解微分方程% f為一階常微分方程的一般表達(dá)式的右端函數(shù)；% *0,y0為初始條件% *n為取值圍的一個(gè)端點(diǎn)；% N為區(qū)間的個(gè)數(shù)；% *為求解微分方程組的值*=zeros(1,N+1); %*為*n構(gòu)成的向量y=zeros(1,N+1); %y為Yn構(gòu)成的向量*(1)=*0;y(1)=y0;h=(*n-*0)/N;for n=1:N *(n+1)=

29、*(n)+h; y(n+1)=y(n)+h*feval(f,*(n),y(n);endT=*,yplot(*,y,b*:);legend(Eule求得的數(shù)值解);在MATLAB命令窗口輸入clear all *=Euler (*,y) y-2*/y,0,1,1,10)運(yùn)行程序，輸出結(jié)果如下：表3.1 Euler法求解一階常微分方程的數(shù)值解*yEuler0 1.00000.1000 1.10000.2000 1.19180.3000 1.27740.4000 1.35820.5000 1.43510.6000 1.50900.7000 1.58030.8000 1.64980.9000 1.71

30、781.0000 1.7848圖3.13.3隱式Euler法的編程實(shí)現(xiàn)隱式Euler法算法見(jiàn)附錄24、用隱式Euler法求解一階常微分方程的數(shù)值解.編寫隱式Euler法的MATLAB程序代碼如下：function T=diEuler(f,*0,y0,*n,N)% diEuler.m為隱式Euler法求微分方程的數(shù)值解% f為一階常微分方程的一般表達(dá)式的右端函數(shù)；% *0,y0為初始條件% *n為取值圍的一個(gè)端點(diǎn)；% n為區(qū)間的個(gè)數(shù)；% *為求解微分方程組的值*=zeros(1,N+1); %*為*n構(gòu)成的向量y=zeros(1,N+1); %y為Yn構(gòu)成的向量*(1)=*0;y(1)=y0;h

31、=(*n-*0)/N;for n=1:N%用迭代法求y(n+1) *(n+1)=*(n)+h; z0=y(n)+h*feval(f,*(n),y(n); for k=1:3 z1=y(n)+h*feval(f,*(n+1),z0); if abs(z1-z0)1e-3; break; end z0=z1; end y(n+1)=z1;endT=*,yplot(*,y,b*:);legend(diEuler求得的數(shù)值解);在MATLAB命令窗口輸入clear all*=diEuler(*,y) y-2*/y,0,1,1,10)運(yùn)行程序，輸出結(jié)果如下：表3.2 隱式Euler法求解一階常微分方程的

32、數(shù)值解* y(diEuler)0 1.0000 0.1000 1.0909 0.2000 1.1743 0.3000 1.2517 0.4000 1.3237 0.5000 1.3910 0.6000 1.4540 0.7000 1.5128 0.8000 1.5676 0.9000 1.6183 1.0000 1.6647圖3.23.4改進(jìn)Euler法的編程實(shí)現(xiàn)改進(jìn)Euler法算法見(jiàn)附錄35、用改進(jìn)Euler法求解一階常微分方程的數(shù)值解編寫改進(jìn)Euler法的MATLAB程序代碼如下：function T=TranEuler(f,*0,y0,*n,N)% TranEuler.m函數(shù)為用改進(jìn)Eu

33、ler法求解微分方程% f為一階常微分方程的一般表達(dá)式的右端函數(shù)；% *0,y0為初始條件% *n為取值圍的一個(gè)端點(diǎn)；% N為區(qū)間的個(gè)數(shù)；% *為求解微分方程組的值*=zeros(1,N+1); %*為*n構(gòu)成的向量y=zeros(1,N+1); %y為Yn構(gòu)成的向量*(1)=*0;y(1)=y0;h=(*n-*0)/N;for n=1:N *(n+1)=*(n)+h; z0=y(n)+h*feval(f,*(n),y(n); y(n+1)=y(n)+h/2*(feval(f,*(n),y(n)+feval(f,*(n+1),z0);endT=*,yplot(*,y,b*:);legend(E

34、uler求得的數(shù)值解);clear all*=TranEuler(*,y) y-2*/y,0,1,1,10)運(yùn)行程序，輸出結(jié)果如下：表3.3 改進(jìn)Euler法求解一階常微分方程的數(shù)值解*y0 1.0000 0.1000 1.0959 0.2000 1.1841 0.3000 1.2662 0.4000 1.3434 0.5000 1.4164 0.6000 1.4860 0.7000 1.5525 0.8000 1.6165 0.9000 1.6782 1.0000 1.7379圖3.3我們已經(jīng)介紹完了關(guān)于微分方程的Euler法、隱式Euler法、改進(jìn)Euler法的MATLAB的編程過(guò)程以及計(jì)

35、算結(jié)果，但是以上三種方法不止一種MATLAB程序，也許有更簡(jiǎn)單的方法在等著讀者去發(fā)現(xiàn)。通過(guò)輸出的數(shù)據(jù)看來(lái)，每種算法得出的結(jié)果各不一樣，則我們終究應(yīng)該選擇哪一種算法更適宜呢？讓我們把以上三種方法的輸出結(jié)果放在一起作為比照，相信你很快會(huì)得你的結(jié)論。首先我們先編寫一個(gè)算法比較程序：clear all%bijiao.m一階微分方程數(shù)值解比照程序%輸入導(dǎo)函數(shù)，*0，y0，*n，NS=(*,y) y-2*/y;R=Euler(S,0,1,1,10);R1=diEuler(S,0,1,1,10);R2=TranEuler(S,0,1,1,10);hold onplot(R(:,1),R(:,2),r*-.)

36、hold onplot(R1(:,1),R1(:,2),g+-)hold onplot(R2(:,1),R2(:,2),b*:)legend(Euler法求得的解,diEuler法求得的解,TranEuler法求得的解)運(yùn)行程序輸出結(jié)果如下：表3.4 比照Euler法，idEuler法以及TranEuler法的數(shù)值解* y(Euler) y(diEuler) y(TranEuler) 0 1.0000 1.0000 1.0000 0.1000 1.1000 1.0909 1.0959 0.2000 1.1918 1.1743 1.1841 0.3000 1.2774 1.2517 1.2662

37、 0.4000 1.3582 1.3237 1.3434 0.5000 1.4351 1.3910 1.4164 0.6000 1.5090 1.4540 1.4860 0.7000 1.5803 1.5128 1.5525 0.8000 1.6498 1.5676 1.6165 0.9000 1.7178 1.6183 1.6782 1.0000 1.7848 1.6647 1.7379圖3.4我們可以直觀的觀察到，在步長(zhǎng)時(shí)，Euler法與改進(jìn)Euler法的曲線貼合較近，隱式Euler法與改進(jìn)Euler法的曲線貼合較遠(yuǎn)。讓我們看看當(dāng)步長(zhǎng)分別為；時(shí)的曲線貼合情況曲線圖像見(jiàn)附錄7。通過(guò)比照可以說(shuō)

38、明，當(dāng)步長(zhǎng)取得較大時(shí)，Euler法的精度不夠高。當(dāng)步長(zhǎng)變小時(shí)，Euler精度會(huì)有所提高。隱式Euler與Euler法情況相似，而改進(jìn)Euler法的精度比前兩種方法的精度都要高，所以當(dāng)時(shí)，隱式Euler與Euler法的數(shù)值解曲線都在向改進(jìn)Euler法的曲線貼合，誤差不斷減??；當(dāng)時(shí)，這三種方法之間的誤差可以忽略不計(jì)。所以我們?cè)诰幊逃?jì)算一階微分方程的時(shí)候，當(dāng)步長(zhǎng)當(dāng)時(shí)，我們應(yīng)該選用改進(jìn)Euler法來(lái)求數(shù)值解；當(dāng)時(shí)，三種方法我們可以任選一種使用。另外當(dāng)步長(zhǎng)取定以后，步數(shù)越多，誤差也會(huì)越大。3.5二階Runge-Kutta法的編程實(shí)現(xiàn)二階Runge-Kutta法算法見(jiàn)附錄46、用二階Runge-Kutta

39、法來(lái)求解一階常微分方程的數(shù)值解編寫二階Runge-Kutta法的MATLAB程序代碼如下：function R=RK2(f,*0,y0,*n,N)% RK2.m函數(shù)為用二階Runge-Kutta法求微分方程的解% f為微分方程% *0,y0為左右端點(diǎn)% *n為給定的初始值% N為給定迭代步長(zhǎng)；% R為求微分方程的解h=(y0-*0)/N;*=zeros(1,N+1);y=zeros(1,N+1);T=*0:h:y0;y(1)=*n;for n=1:N k1=h*feval(f,T(n),y(n); k2=h*feval(f,T(n)+2*h/3,y(n)+2*h*k1/3); y(n+1)=y

40、(n)+h*(k1+3*k2)/4;endR=T,y編寫二階Heun法函數(shù)圖像的MATLAB程序代碼如下：S=(*,y) *-y+1;t1,y1=ode23(S,0:0.1:1,1);R=RK2(S,0,1,1,10)plot(R(:,1),R(:,2),r)hold onplot(t1,y1,b* -.)legend(二階Heun方法求得的解,ode23求得的解)hold offK=t1,y1將這兩個(gè)M文件分別命名為RK2.m和RK2_Heun.m放在一個(gè)文件夾中，然后直接運(yùn)行RK2_Heun.m程序，輸出結(jié)果如下：表3.5二階Runge-Kutta法求解一階常微分方程的數(shù)值解* y(RK2

41、) y(ode23)0 1.0000 1.00000.1000 1.0000 1.00030.2000 1.0003 1.00250.3000 1.0009 1.00840.4000 1.0021 1.01940.5000 1.0041 1.03690.6000 1.0071 1.06240.7000 1.0112 1.09680.8000 1.0167 1.14130.9000 1.0238 1.19681.0000 1.0325 1.2642圖3.53.6三階Runge-Kutta法的編程實(shí)現(xiàn)三階Runge-Kutta法算法見(jiàn)附錄57、用三階Runge-Kutta法來(lái)求解一階常微分方程的數(shù)

42、值解編寫三階Runge-Kutta法的MATLAB程序代碼如下：function R=RK3(f,*0,y0,*n,N)% RK3.m函數(shù)為用三階Runge-Kutta法求微分方程的解% f為微分方程% *0,y0為左右端點(diǎn)% *n為給定的初始值% N為給定迭代步長(zhǎng)；% R為求微分方程的解h=(y0-*0)/N;*=zeros(1,N+1);y=zeros(1,N+1);T=*0:h:y0;y(1)=*n;for n=1:N k1=h*feval(f,T(n),y(n); k2=h*feval(f,T(n)+h/2,y(n)+k1/2); k3=h*feval(f,T(n)+3*h/4,y(n

43、)+3*h/4*k2); y(n+1)=y(n)+(2*k1+3*k2+4*k3)/9;endR=T,y編寫三階Runge-Kutta法函數(shù)圖像的MATLAB程序代碼如下：S=(*,y) *2-y+1;t1,y1=ode23(S,0:0.1:1,1)R=RK3(S,0,1,1,10)plot(R(:,1),R(:,2),r)hold onplot(t1,y1,b*-.)legend(三階RK法求得的解,ode23求得的解)hold off將這兩個(gè)M文件分別命名為RK3.m和RK3_1.m放在一個(gè)文件夾中，然后直接運(yùn)行RK3_1.m程序，輸出結(jié)果如下：表3.6 三階Runge-Kutta法求解一

44、階常微分方程的數(shù)值解*yRK3 yode230.1000 1.0003 1.00030.4000 1.0199 1.01940.6000 1.0641 1.06240.7000 1.0994 1.09680.8000 1.1450 1.14130.9000 1.2018 1.19681.0000 1.2707 1.2642圖3.63.7四階Runge-Kutta法的編程實(shí)現(xiàn)四階Runge-Kutta法算法見(jiàn)附錄68、用四階Runge-Kutta法來(lái)求解一階常微分方程的數(shù)值解編寫四階Runge-Kutta法的MATLAB程序代碼如下：function R=RK4(f,*0,y0,*n,N)% R

45、K4.m函數(shù)為用四階Runge-Kutta法求微分方程的解% f為微分方程% *0,y0為左右端點(diǎn)% *n為給定的初始值% N為給定迭代步長(zhǎng)；% R為求微分方程的解h=(y0-*0)/N;*=zeros(1,N+1);y=zeros(1,N+1);T=*0:h:y0;y(1)=*n;for n=1:N k1=h*feval(f,T(n),y(n); k2=h*feval(f,T(n)+h/2,y(n)+k1/2); k3=h*feval(f,T(n)+h/2,y(n)+k2/2); k4=h*feval(f,T(n)+h,y(n)+k3); y(n+1)=y(n)+(k1+2*k2+2*k3+

46、k4)/6;endR=T,y編寫四階Runge-Kutta法函數(shù)圖像的MATLAB程序代碼如下：S=(*,y) *2-y+1;t1,y1=ode45(S,0:0.1:1,1)R=RK4(S,0,1,1,10)plot(R(:,1),R(:,2),r)hold onplot(t1,y1,b*-.)legend(四階RK法求得的解,ode45求得的解)hold off將這兩個(gè)M文件分別命名為RK4.m和RK4_1.m放在一個(gè)文件夾中，然后直接運(yùn)行RK4_1.m程序，輸出結(jié)果如下：表3.7 四階Runge-Kutta法求解一階微分方程的數(shù)值解* y(RK4) y(ode45) 0 1.0000 1.

47、0000 0.1000 1.0003 1.0003 0.2000 1.0025 1.0025 0.3000 1.0084 1.0084 0.4000 1.0194 1.0194 0.5000 1.0369 1.0369 0.6000 1.0624 1.0624 0.7000 1.0968 1.0968 0.8000 1.1413 1.1413 0.9000 1.1969 1.1969 1.0000 1.2642 1.2642圖 3.7我們可以看到，用四階Runge-Kutta法求一階微分方程的數(shù)值解曲線和ode45求得的數(shù)值解曲線完全重合，說(shuō)明在給定步長(zhǎng)的情況下，四階Runge-Kutta法階

48、的精度到達(dá)了我們求一階微分方程的精度要求，與準(zhǔn)確解之間的誤差非常小。接下來(lái)我們可把以上三種Runge-Kutta方法進(jìn)展比較一下，看看對(duì)于同樣的微分方程各數(shù)值方法的數(shù)值解之間的誤差是多少，以幫助我們更好更恰當(dāng)?shù)倪x擇數(shù)值方法。編寫一階微分方程數(shù)值解比照的MATLAB程序代碼如下：clear all%duibi.m一階微分方程數(shù)值解比照程序%輸入導(dǎo)函數(shù)，*0，y0，*n，h，NS=(*,y) *-y+1;t1,y1=ode23(S,0:0.1:1,1);t2,y2=ode45(S,0:0.1:1,1)R=RK2(S,0,1,1,10);K=RK3(S,0,1,1,10);G=RK4(S,0,1,1

49、,10);plot(R(:,1),R(:,2),r,K(:,1),K(:,2),c,G(:,1),G(:,2),r)hold onplot(t1,y1,b* -,t2,y2,m+ :)legend(二階Heun方法求得的解,三階RK法求得的解,四階RK法求得的解,ode23求得的解,ode45求得的解)hold off將程序命名為duibi.m，和RK2.m，RK3.m，RK4.m放在一起，運(yùn)行duibi.m程序，輸出函數(shù)圖像結(jié)果請(qǐng)看圖 3.7。表3.8一階微分方程五種數(shù)值解比照表* y(RK2) y(RK3) y(RK4) y(ode23) y(ode45)0 1.0000 1.00001.

50、0000 1.0000 1.00000.1000 1.00001.00031.0003 1.0003 1.00030.2000 1.00031.00261.0025 1.0025 1.00250.3000 1.00091.0086 1.0084 1.0084 1.00840.4000 1.0021 1.0199 1.0194 1.0194 1.01940.5000 1.00411.0380 1.0369 1.0369 1.03690.6000 1.00711.0641 1.0624 1.0624 1.06240.7000 1.01121.0994 1.0968 1.0968 1.09680.8

51、000 1.01671.1450 1.1413 1.1413 1.14130.9000 1.0238 1.2018 1.1969 1.1968 1.19691.0000 1.0325 1.2707 1.2642 1.2642 1.2642圖3.8我們通過(guò)觀察圖 3.7我們可以發(fā)現(xiàn)二階Heun法的數(shù)值解曲線與其它數(shù)值解曲線偏離程度較大，而且是增大趨勢(shì)，而其它數(shù)值解曲線相互之間偏離程度小，其中四階Runge-Kutta法與ode45法的數(shù)值解曲線重合，與ode23法的數(shù)值解曲線幾乎重合。在MATLAB中，函數(shù)ode23和ode45采用的是變步長(zhǎng)的Runge-Kutta法，并且ode45的精度更高一

52、些。ode23適合的常微分方程問(wèn)題類型為低精度容差或適度非剛性問(wèn)題；ode45適合的常微分方程問(wèn)題類型為非剛性方程。所以在精度誤差的允許圍，我們優(yōu)先選擇四階Runge-Kutta法或者函數(shù)ode45來(lái)解決一階微分方程的初值問(wèn)題。則在編程實(shí)現(xiàn)的時(shí)候，我們?cè)撊绾芜x擇適當(dāng)?shù)腞unge-Kutta法呢？我們同樣通過(guò)改變步長(zhǎng)的方法來(lái)研究該問(wèn)題，分別取步長(zhǎng)為；，觀察曲線的貼合情況Runge-Kutta法曲線比照見(jiàn)附錄8。當(dāng)步長(zhǎng)在變小的過(guò)程中，RK3與RK4的曲線越來(lái)越貼近，重合程度越來(lái)越大，并且與RK2偏離程度也越來(lái)越大，說(shuō)明RK2的精度比較低的，當(dāng)步長(zhǎng)之后，RK2的精度可能會(huì)比改進(jìn)的Euler法還要低，

53、所以我們?cè)诰幊虒?shí)現(xiàn)一階微分方程的數(shù)值解并對(duì)步長(zhǎng)嚴(yán)格控制的時(shí)候，我們應(yīng)該優(yōu)先選擇三階或者四階的Runge-Kutta法，控制數(shù)值解與準(zhǔn)確解之間的誤差。結(jié) 論本文探究的一階微分方程的MATLAB數(shù)值解法有Euler法、隱式Euler、改進(jìn)Euler法、二階Runge-Kutta法、三階Runge-Kutta法以及四階Runge-Kutta法。采用比照分析的方法探究了以上方法的精度問(wèn)題，如何編程實(shí)現(xiàn)算法，以及如何選擇程序求解一階微分方程初值問(wèn)題。在給定步長(zhǎng)的微分方程初值問(wèn)題中，得出以下結(jié)論：利用Euler法可以求解簡(jiǎn)單的微分方程初值問(wèn)題，Euler法是最簡(jiǎn)單的單步法，MATLAB編程實(shí)現(xiàn)操作簡(jiǎn)單，很

54、好的表達(dá)了算法的決策思想，在步長(zhǎng)充分小的時(shí)候，Euler法才可用，一般不用于實(shí)際計(jì)算，但是使用方便，是一階精度的算法；隱式Euler法也被稱為后退的Euler法，隱式Euler法的算法在穩(wěn)定性方面表現(xiàn)出色，但是迭代次數(shù)比Euler法多，計(jì)算過(guò)程比較久，在特殊場(chǎng)合人們也會(huì)使用該算法來(lái)計(jì)算；改進(jìn)Euler法相比前面的兩種算法來(lái)說(shuō)，算法的精度是二階的，具有良好的穩(wěn)定性，受步長(zhǎng)的影響不大，被稱為預(yù)測(cè)矯正系統(tǒng)，MATLAB編程實(shí)現(xiàn)過(guò)程簡(jiǎn)單，人們通常會(huì)選用該算法來(lái)求解一階微分方程的初值問(wèn)題；二階Runge-Kutta法的算法精度是二階的，二階Runge-Kutta法可以推導(dǎo)出改進(jìn)的Euler法，也可以換算

55、成文中有介紹的中點(diǎn)方法和Heun方法；三階Runge-Kutta法的算法精度是三階的，穩(wěn)定性也比教好；四階Runge-Kutta法的算法精度是四階的，是以上的算法中階的精度最高的算法，穩(wěn)定性也以上的算法中最好的，四階Runge-Kutta法的公式比較多，人們經(jīng)常選擇經(jīng)典的四階Runge-Kutta法來(lái)解決一階微分方程的初值問(wèn)題。Runge-Kutta法的編程實(shí)現(xiàn)復(fù)雜程度比Euler法的要難一些，需要考慮Runge-Kutta算法的步長(zhǎng)的取值，并且由于迭代的次數(shù)較多，所以算法的計(jì)算量大；單步法的收斂性與穩(wěn)定性：針對(duì)給定步長(zhǎng)的一階微分方程初值問(wèn)題，三種Euler法和二到四階Runge-Kutta法

56、都具有收斂性；在實(shí)際應(yīng)用中，選擇數(shù)值方法的步長(zhǎng)時(shí)，我們應(yīng)該考慮穩(wěn)定性，步長(zhǎng)的取值變化應(yīng)該在絕對(duì)穩(wěn)定區(qū)間，以保證數(shù)值解與準(zhǔn)確解之間的誤差。參 考 文 獻(xiàn)1王高雄，朱四銘，周之銘，王壽松. 常微分方程M. ：高等教育，20062C.W.吉爾. 常微分方程初值問(wèn)題的數(shù)值解法M. 費(fèi)景高 譯. ：科學(xué)，19783德豐. MATLAB數(shù)值分析M. 第二版.：機(jī)械工業(yè)，20124樓順天，于衛(wèi)，閆華梁. MATLAB程序設(shè)計(jì)語(yǔ)言M. ：電子科技大學(xué)，19975夏省祥，于正文. 常用數(shù)值算法及其MATLAB實(shí)現(xiàn)M. ：清華大學(xué)，20146白峰杉. 數(shù)值計(jì)算引論M.第二版. ：高等教育，20107白峰杉，慶揚(yáng)，關(guān)治. 數(shù)值計(jì)算原理M.：清華大學(xué)，20008文林. 數(shù)學(xué)史概論M. ：高等教育，20029杜瑞芝. 數(shù)學(xué)史辭典M. ：教育，200010VictorJ.katz. 數(shù)學(xué)史通論M. 文林 譯. ：高等教育，200411Gerald Recktenwald. 數(shù)值方法和MATLAB實(shí)現(xiàn)與應(yīng)用M.伍衛(wèi)國(guó)，萬(wàn)群，輝譯.：機(jī)械工業(yè)，200412吳開(kāi)騰，譚燕梅，