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

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

2、體的微分方程來研究以上算法的編程實現(xiàn)，通過MATBLAB求解具體的一階微分方程來探究以上方法的優(yōu)缺點，通過圖表來分析得出結論：改進Euler法和四階Runge-Kutta法的階的精度較高，具有較好的算法穩(wě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常微分方程的開展簡史 PAGEREF _Toc21991 21.3 國外研究現(xiàn)狀 PAGEREF _Toc15673 21.4 研究主要容及研究方案 PAGEREF _Toc11775 31.4.1 研究的主要容 PAGEREF _Toc1002 3 研究方案 PAGEREF _Toc297

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

27、的末尾定義該方程中的自變量為“*，否則系統(tǒng)會把“t認為自變量，從而造成錯誤，所以要引起注意。我們通過兩個例子來了解一下dsolve函數(shù)的使用方法：1、計算微分方程的解在MATLAB命令窗口輸入clear allz=dsolve(Dy=-*,D*=y)運行程序輸出結果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、計算微分方程的通解在命令窗口輸入以下命令clear ally=dsolve(Dy=*2*y,*)輸出結果y = C2*e*p(*2)所

28、以，通解為3.2Euler法的編程實現(xiàn)Euler法算法見附錄13、用Euler法求解一階常微分方程的數(shù)值解編寫Euler法的MATLAB程序代碼如下：function T= Euler(f,*0,y0,*n,N)% Euler.m函數(shù)為用Euler法求解微分方程% f為一階常微分方程的一般表達式的右端函數(shù)；% *0,y0為初始條件% *n為取值圍的一個端點；% N為區(qū)間的個數(shù)；% *為求解微分方程組的值*=zeros(1,N+1); %*為*n構成的向量y=zeros(1,N+1); %y為Yn構成的向量*(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)運行程序，輸出結果如下：表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法的編程實現(xiàn)隱式Euler法算法見附錄24、用隱式Euler法求解一階常微分方程的數(shù)值解.編寫隱式Euler法的MATLAB程序代碼如下：function T=diEuler(f,*0,y0,*n,N)% diEuler.m為隱式Euler法求微分方程的數(shù)值解% f為一階常微分方程的一般表達式的右端函數(shù)；% *0,y0為初始條件% *n為取值圍的一個端點；% n為區(qū)間的個數(shù)；% *為求解微分方程組的值*=zeros(1,N+1); %*為*n構成的向量y=zeros(1,N+1); %y為Yn構成的向量*(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)運行程序，輸出結果如下：表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改進Euler法的編程實現(xiàn)改進Euler法算法見附錄35、用改進Euler法求解一階常微分方程的數(shù)值解編寫改進Euler法的MATLAB程序代碼如下：function T=TranEuler(f,*0,y0,*n,N)% TranEuler.m函數(shù)為用改進Eu

33、ler法求解微分方程% f為一階常微分方程的一般表達式的右端函數(shù)；% *0,y0為初始條件% *n為取值圍的一個端點；% N為區(qū)間的個數(shù)；% *為求解微分方程組的值*=zeros(1,N+1); %*為*n構成的向量y=zeros(1,N+1); %y為Yn構成的向量*(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)運行程序，輸出結果如下：表3.3 改進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)介紹完了關于微分方程的Euler法、隱式Euler法、改進Euler法的MATLAB的編程過程以及計

35、算結果，但是以上三種方法不止一種MATLAB程序，也許有更簡單的方法在等著讀者去發(fā)現(xiàn)。通過輸出的數(shù)據(jù)看來，每種算法得出的結果各不一樣，則我們終究應該選擇哪一種算法更適宜呢？讓我們把以上三種方法的輸出結果放在一起作為比照，相信你很快會得你的結論。首先我們先編寫一個算法比較程序：clear all%bijiao.m一階微分方程數(shù)值解比照程序%輸入導函數(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法求得的解)運行程序輸出結果如下：表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我們可以直觀的觀察到，在步長時，Euler法與改進Euler法的曲線貼合較近，隱式Euler法與改進Euler法的曲線貼合較遠。讓我們看看當步長分別為；時的曲線貼合情況曲線圖像見附錄7。通過比照可以說

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

39、法來求解一階常微分方程的數(shù)值解編寫二階Runge-Kutta法的MATLAB程序代碼如下：function R=RK2(f,*0,y0,*n,N)% RK2.m函數(shù)為用二階Runge-Kutta法求微分方程的解% f為微分方程% *0,y0為左右端點% *n為給定的初始值% N為給定迭代步長；% 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將這兩個M文件分別命名為RK2.m和RK2_Heun.m放在一個文件夾中，然后直接運行RK2_Heun.m程序，輸出結果如下：表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法的編程實現(xiàn)三階Runge-Kutta法算法見附錄57、用三階Runge-Kutta法來求解一階常微分方程的數(shù)

42、值解編寫三階Runge-Kutta法的MATLAB程序代碼如下：function R=RK3(f,*0,y0,*n,N)% RK3.m函數(shù)為用三階Runge-Kutta法求微分方程的解% f為微分方程% *0,y0為左右端點% *n為給定的初始值% N為給定迭代步長；% 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將這兩個M文件分別命名為RK3.m和RK3_1.m放在一個文件夾中，然后直接運行RK3_1.m程序，輸出結果如下：表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法的編程實現(xiàn)四階Runge-Kutta法算法見附錄68、用四階Runge-Kutta法來求解一階常微分方程的數(shù)值解編寫四階Runge-Kutta法的MATLAB程序代碼如下：function R=RK4(f,*0,y0,*n,N)% R

45、K4.m函數(shù)為用四階Runge-Kutta法求微分方程的解% f為微分方程% *0,y0為左右端點% *n為給定的初始值% N為給定迭代步長；% 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將這兩個M文件分別命名為RK4.m和RK4_1.m放在一個文件夾中，然后直接運行RK4_1.m程序，輸出結果如下：表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ù)值解曲線完全重合，說明在給定步長的情況下，四階Runge-Kutta法階

48、的精度到達了我們求一階微分方程的精度要求，與準確解之間的誤差非常小。接下來我們可把以上三種Runge-Kutta方法進展比較一下，看看對于同樣的微分方程各數(shù)值方法的數(shù)值解之間的誤差是多少，以幫助我們更好更恰當?shù)倪x擇數(shù)值方法。編寫一階微分方程數(shù)值解比照的MATLAB程序代碼如下：clear all%duibi.m一階微分方程數(shù)值解比照程序%輸入導函數(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放在一起，運行duibi.m程序，輸出函數(shù)圖像結果請看圖 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我們通過觀察圖 3.7我們可以發(fā)現(xiàn)二階Heun法的數(shù)值解曲線與其它數(shù)值解曲線偏離程度較大，而且是增大趨勢，而其它數(shù)值解曲線相互之間偏離程度小，其中四階Runge-Kutta法與ode45法的數(shù)值解曲線重合，與ode23法的數(shù)值解曲線幾乎重合。在MATLAB中，函數(shù)ode23和ode45采用的是變步長的Runge-Kutta法，并且ode45的精度更高一

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

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

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

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

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