常微分方程常用數(shù)值解法

1、第一章 緒論1.1 引言常微分方程是現(xiàn)代數(shù)學(xué)的一個重要分支，是人們解決各種實際問題的有效工具。微分方程的理論和方法從17世紀末開始發(fā)展起來，很快成了研究自然現(xiàn)象的強有力工具，在17到18世紀，在力學(xué)、天文、科學(xué)技術(shù)、物理中，就已借助微分方程取得了巨大的成就。1864年Leverrer根據(jù)這個方程預(yù)見了海王星的存在，并確定出海王星在天空中的位置?，F(xiàn)在，常微分方程在許多方面獲得了日新月異的應(yīng)用。這些應(yīng)用也為常微分方程的進一步發(fā)展提供了新的問題，促使人們對微分方程進行更深入的研究，以便適應(yīng)科學(xué)技術(shù)飛速發(fā)展的需要。研究常微分方程常用數(shù)值解是數(shù)學(xué)工作者的一項基本的且重要的工作。在國內(nèi)外眾多數(shù)學(xué)家的不懈努

2、力，使此學(xué)科基本上形成了一套完美的體系。微分方程的首要問題是如何求一個給定方程的通解或特解。到目前為止，人們已經(jīng)對許多微分方程得出了求解的一般方法。由于在生產(chǎn)實際和科學(xué)研究中所遇到的微分方程問題比較復(fù)雜，使這些問題的解即使能求出解析表達式，也往往因計算量太大而難于求出，而對于一些典型的微分方程則可以運用基本方法求出其解析解，并可以根據(jù)初值問題的條件把其中的任意常數(shù)確定下來。由于求通解存在許多困難，人們就開始研究帶某種定解條件的特解。首先是Cauchy對微分方程初始解的存在惟一性進行了研究。目前解的存在惟一性、延拓性、大范圍的存在性以及解對初始解和參數(shù)的延續(xù)性和可微性等理論問題都已發(fā)展成熟。與此

3、同時，人們開始采取各種近似方法來求微分方程的特解，例如求微分方程數(shù)值解的Euler折線法、Runge-Kutta法等，可以求得若干個點上微分方程的近似解。最后，由于當代高科技的發(fā)展為數(shù)學(xué)的廣泛應(yīng)用和深入研究提供了更好的手段。用計算機結(jié)合Matlab軟件求方程的精確解、近似解，對解的性態(tài)進行圖示和定性、穩(wěn)定性研究都十分方便有效。本章先介紹常微分的一般概念、導(dǎo)出微分方程的一些典型例子及求解微分方程的思路分析。從而得到常微分方程的常用數(shù)值解法。1.2 常微分方程的概念1.常微分方程的定義含有未知量的等式稱為方程，它表達了未知量所必須滿足的某些條件。一般說來，凡含有自變量、未知函數(shù)以及未知函數(shù)的導(dǎo)數(shù)或

4、微分的方程稱為微分方程。如果微分方程中的未知函數(shù)只依賴于一個自變量，則稱為常微分方程；如果未知函數(shù)依賴于兩個或多個的自變量，并且在方程中出現(xiàn)偏導(dǎo)數(shù)，則稱為偏微分方程。在在一個微分方程中，未知函數(shù)最高階導(dǎo)數(shù)的階數(shù)，稱為方程的階數(shù)。如果一個微分方程關(guān)于未知函數(shù)及其各階導(dǎo)數(shù)都是線性的，則稱它為線性微分方程，否則稱之為非線性微分方程。本論文主要介紹常微分方程，也簡稱微分方程。 以為未知函數(shù)，為自變量的一階常微分方程的一般形式可表示為微分方程， (1.2.1) 將（1.2.1）中解出，則得到方程 (1.2.2)或 (1.2.3)也稱（1.2.1）為一階隱式微分方程，（1.2.2）為一階顯式微分方程，（1

5、.2.3）為一階微分方程的微分形式。 階隱式方程的一般形式為 (1.2.4) 階顯式方程的一般形式 (1.2.5)方程（1.2.4）中，如果函數(shù)對未知函數(shù)和它的各階導(dǎo)數(shù)都是一次的，則稱其為線性常微分方程，否則，稱其為非線性微分方程。以為未知函數(shù)，為自變量的階線性微分方程具有如下形式： (1.2.6) 2.常微分方程的解 設(shè)函數(shù)是定義在區(qū)間上的階可微導(dǎo)數(shù)。如果把代入方程（1.2.4）后能使其成為恒等式，即 則稱是微分方程（1.2.4）在區(qū)間上的一個解。例如，是微分方程在的一個解.是微分方程在區(qū)間的一個解.如果關(guān)系式?jīng)Q定的隱函數(shù)是方程（1.2.4）的解，則我們稱是（1.2.4）的隱式解。例如一階微

6、分方程 有隱式解 我們把含有個相互獨立的任意常數(shù)的解 稱為階微分方程（1.2.4）的通解。在通解之中當一組任意常數(shù)確定時，所得到確定的解稱為特解。例如，是二階線性方程的通解，而都是其特解，其中是任意常數(shù)。一般地，方程的特解可由其通解中任意常數(shù)取確定的常數(shù)導(dǎo)出，且方程的通解不一定表示方程的所有解。3. 常微分方程初值問題為了確定微分方程的一個特解，我們可以給出這個微分方程的所滿足的定解條件，常見的定解條件是初始條件，即方程（1.2.4）在某一點所滿足的條件： （1.2.7）微分方程（1.2.4）連同初始條件（1.2.7）一起稱為初始值問題。 一階常微分方程初值問題是求解函數(shù)，且滿足 （1.2.8

7、）其中為已知函數(shù)，為給定的初值。定理1.2.1 假設(shè)在矩形區(qū)域內(nèi)連續(xù)，且關(guān)于變元Lipschitz連續(xù)，即存在正常數(shù)，使得對任意及，成立不等式 (1.2.9)其中常數(shù)稱為Lipschitz常數(shù),則處置問題存在惟一解，連續(xù)依賴于初值由常微分方程的基本理論，我們有：換句話說，若是如下初值問題 (1.2.10)的解，則 因此問題（1.2.8）是適定的。定理1.2.2 當假定函數(shù)滿足定理1.2.1中的條件時，也即微分方程（1.2.8）的解是適定的。 綜上，高階常微分方程初值問題一般為 (1.2.11)其中為給定值。引進新的變量函數(shù) (1.2.12)則初值問題（1.2.11）化成了一階常微分方程組初值問

8、題 (1.2.13)通過求解（1.2.13）得到（1.2.11）的解1.3 常微分方程常用數(shù)值解法的思路一般說來，對于某些典型微分方程可以通過初等積分法求解，而對于復(fù)雜微分方程，需要得到解在若干個點上的近似值或者便于計算的近似表達式即可。本論文講研究常微分方程常用數(shù)值解法及其matlab程序設(shè)計。對于一階常微分方程的初值問題為 (1.3.1)求微分方程初值問題（1.3.1）的數(shù)值解，就是求解函在一系列離散點上精確值的近似值。由于計算的復(fù)雜及繁冗，在計算常微分方程初值問題的數(shù)值解是一般使用Matlab軟件進行編程計算，按如下步驟：（1）引入點列，其中，稱為步長。為了便于使用計算機進行編程計算，一

9、般取步長為定值，即， (2.2)（2）使用常微分方程常用數(shù)值解方法，即輸入由計算出的遞推公式。（3）利用（2）中的公式逐步求出近似解第二章 單步法所謂單步法是指這類方法在計算時，只用到前一步的值，然后逐步往下計算。單步法的精度不高，這是它的缺點。但我們只要給定初值就可以計算，所以還是比較方便的。以下我們介紹兩種典型的單步法：Euler法和Runge-Kutta法。2.1 Euler法為了計算出初始值問題 （2.1.1） 在這一區(qū)間上若干個點上解的近似解，我們將此區(qū)間等分，令為（2.1.1）的解，即當時精確值的近似解。我們的任務(wù)是要用盡可能簡單的方法求出與盡可能接近的值。Euler法又稱Eule

10、r折線法，Euler折線法是數(shù)值解法最簡單的一種。Euler折線法的基本思想是利用微分中值定理對（2.1.1）的解函數(shù)進行近似。由于是可微函數(shù)，故 其中，是介于和之間的一個值。當比較小時，和相差不大，故用代替上式中的，就得到了的近似值 當知道以后，再用類似方法來求在點的近似值 （2.1.2）這樣就得到了（2.1.1）的解在點的近似解，從幾何上看，Euler折線法就是在局部范圍內(nèi)用切線上的值去替代解曲線上的值。Euler公式在實際計算時較少采用，但由于它的結(jié)構(gòu)簡單，易于分析，在理論上具有非常重要的意義。2.1.1 向前Euler法Euler折線法的計算量小，但誤差大，后來人們對Euler折線法進

11、行了改進：由方程（2.1.1），在節(jié)點處成立 （2.1.3）將在處利用Taylor展開式展開，得 其中將將式（2.1.3）代入上式得 （2.1.4） 在式（2.1.4）中略去高階項 （2.1.5）并用和分別代替式（2.1.4）中的和，可得Euler公式（也稱向前Euler公式： （2.1.6）2.1.2 向后Euler法向后Euler法和向前Euler法相差不大，只是將在處利用Taylor展開式展開，得 其中將將式（2.1.3）代入上式得 （2.1.7） 在式（2.1.7）中略去高階項，并用和分別代替式（2.1.7）中的和，可得向后Euler公式： （2.1.8）其中，局部截斷誤差為 因此，向

12、后Euler公式同向前Euler公式，為一階方法。 向前Euler公式（2.1.6）和向后Euler公式（2.1.8）相比，向前Euler公式關(guān)于是顯式格式。而向后Euler公式中，隱含在方程中，稱這樣的格式為隱式格式。2.1.3 中點差分公式法若將展開式（2.1.4）和（2.1.7）多展開一項，得 其中將上述兩式相減并整理得 其中略去高階項，并用和代替上式中 和，即得中心差分公式 （2.1.9）其局部截斷誤差為 即中心差分公式為二階方法。2.1.4 梯形公式如果將（2.1.1）中微分方程兩端從到積分，則得 b 如果用梯形求積公式計算上式中的積分，則有 略去高階項可得梯形公式： (2.1.10

13、)其中，局部截斷誤差為 顯然，梯形公式為二階隱式公式。2.1.5 改進Euler法一般來說，隱式格式比顯示格式具有較好的數(shù)值穩(wěn)定性。但是，由于在隱式格式里，往往滿足的是一個非線性方程，因而需要使用迭代法得到的一個近似值，然后代入隱式格式作校正，并以這個校正值作為的近似值。印版地，可以選擇適當?shù)娘@示格式計算預(yù)測值，我們把這樣的格式稱為預(yù)估一校正公式。預(yù)估一校正公式往往既便于計算又具有較好的數(shù)值穩(wěn)定性。如，對于梯形公式，可用向前Euler公式計算出的預(yù)測值，再用梯形公式修正得到： 上述公式為改進Euler公式。 其總體截斷誤差為在步長h取得足夠小的時侯，我們用改進Euler方法解得數(shù)值與理論值非常

14、接近。2.2 Runge-Kutta法Runge-Kutta法是一種間接使用Taylor展開式來獲取高階方法的數(shù)值方法，也是一種特殊的單步法，具有相當實用價值。其基本做法是：首先，用函數(shù)在附近的個點上的函數(shù)值的線性組合來代替的導(dǎo)數(shù)。然后，通過Taylor展開式，確定其中的組合系數(shù)。這樣既可以避免計算高階導(dǎo)數(shù)，又可以獲得較高的方法階。設(shè)為一個正整數(shù)（代表使用的函數(shù)值的個數(shù)）， 和 為一些特定參數(shù)，則數(shù)值方法 (2.2.1)稱為級Runge-Kutta公式。2.2.1 二階Runge-Kutta法顯然，公式（2.2.1）是顯式的。如果對于該代數(shù)方程組的一組解，對應(yīng)的方法階為，則得到一個級階顯式Ru

15、nge-Kutta公式。當時，公式（2.2.1）變?yōu)?其局部截斷誤差為 經(jīng)計算得參數(shù)方程組： （2.2.2） 1.若取則對應(yīng)兩二級二階方法（又稱變形Euler公式）： 2. 若取則對應(yīng)兩級二階方法（又稱兩級Huen公式）： 2.2.2 三階Runge-Kutta法1. 三級Kutta公式（三級三階格式） 2. 三級Huen公式（三級三階格式） 這兩個常用方法在實際使用中能達到較低的精確要求。2.2.3 古典Runge-Kutta法古典Runge-Kutta公式是精度高、最實用的一種常微分方程常用數(shù)值解法。 2.2.4 一級隱式中點公式設(shè)為一個正整數(shù)（代表使用的函數(shù)值的個數(shù)）， 和 為一些特定參

16、數(shù)，則數(shù)值方法 (2.2.2) 稱為級隱式Runge-Kutta公式。當時，有一級隱式中點公式為 2.2.5 二級隱式中點公式當時，有二級隱式中點公式為 相對顯式公式而言，對于同樣的，隱式格式不僅可獲取較高的精度，而且具有較好的數(shù)值穩(wěn)定性。第三章 多步法 所謂多步法指：這類方法在計算時，除了用到前一步的值，之外，還要用到 ，這前面步的值，這個算法的代表就是阿達姆斯(Adams)方法。 本章介紹線性多步法。線性多步法多使用于求解步驟較多的情況。一般地，步線性多步法公式可寫成 （3.1）當時，式（3.1）為顯式格式。否則為隱式格式。1. 顯式格式首先，將常微分方程（2.1.1）兩端在區(qū)間上求定積分

17、，得 如果在上面等式中的定積分用關(guān)于結(jié)點的數(shù)值積分近似，可得如下顯式格式 其中 而為函數(shù)關(guān)于結(jié)點的Lagrange插值函數(shù)的基函數(shù)。 當時，四步Adams顯式格式 其局部截斷誤差為 2.隱式格式如果在等式（3.1）中的定積分用關(guān)于結(jié)點的數(shù)值積分近似，可得如下隱式格式 其中 而為函數(shù)關(guān)于結(jié)點的Lagrange插值函數(shù)的基函數(shù)。當時，三步Adams隱式格式 其局部截斷誤差為 Adams顯式與隱式法相比較而言：同階數(shù)情況下，隱式格式的局部截斷誤差系數(shù)的絕對值比顯式格式的小，隱式格式的穩(wěn)定性相對顯式格式的要大。 總結(jié) 本論文總結(jié)了單步法（向前Euler法，向后Euler法，中心差分法，梯形法，改進Eu

18、ler法，Runge-Kutta法）及線性多步法（Adams法）幾種比較常用的常微分方程數(shù)值解法，并通過數(shù)值算例對同一問題的解的誤差計算比較，總結(jié)各種算法的優(yōu)缺點。一般而言，對Euler法而言，求解常微分方程初值問題數(shù)值時計算簡單、程序簡潔并且計算量小，但是，計算精度卻很低。中心差分法及梯形法較好地提高計算精確度，同時又保持計算簡單，程序簡潔的優(yōu)點。改進的Euler法比Euler法得到的計算結(jié)果要好很多。我們用改進Euler方法解得數(shù)值與理論值非常接近。對Runge-Kutta法來說，越高階的方法計算精度越高，但相應(yīng)的計算量也越大。古典Runge-Kutta公式是精度高、最實用的一種常微分方程

19、常用數(shù)值解法。而線性多步法多使用于求解步驟較多的情況。 我們要把常微分方程常用數(shù)值解法掌握并靈活運用并加以推廣應(yīng)用，來解決我們的實際生活中所遇到問題，為我國的發(fā)展服務(wù)，充分體現(xiàn)了常微分方程常用數(shù)值解法的重大意義。作為一名數(shù)學(xué)工作者，我們有責任也有義務(wù)繼續(xù)探索常微分方程組的數(shù)值解法，并且在不斷探索其解法的前提下，開發(fā)出更好更優(yōu)的軟件來解決實際算法問題。但是，其現(xiàn)有理論也還遠遠不能滿足需要，有待于進一步研究，使這門學(xué)科的理論更加完善。 數(shù)值實例1.分別用改進Euler公式和古典Runge-Kutta公式解常微分方程初值問題取步長，計算結(jié)果保存于下表1中，精確解為。解：fun函數(shù)： function

20、 f=fun(x,y); f=y-2*x/y;精確解： x=0:0.1:1; y=sqrt(1+2*x)y = Columns 1 through 9 1.0000 1.0954 1.1832 1.2649 1.3416 1.4142 1.4832 1.5492 1.6125 Columns 10 through 11 1.6733 1.73211. 應(yīng)用改進Euler法，運用matlab軟件，程序為：format long h=0.1;x=0:h:1.0;y=zeros(size(x);y(1)=1;y1=sqrt(1+2*x);N=length(x);for n=1:N-1 K1=feva

21、l(fun,x(n),y(n); K2=feval(fun,x(n)+0.5*h,y(n)+0.5*h*K1); y(n+1)=y(n)+h*K1;endy2. 應(yīng)用古典Runge-Kutta法，運用matlab軟件，程序為：format longh=0.1;x=0:h:1.0;y=zeros(size(x);y(1)=1;y1=sqrt(1+2*x);N=length(x);for n=1:N-1 K1=feval(fun,x(n),y(n); K2=feval(fun,x(n)+0.5*h,y(n)+0.5*h*K1); K3=feval(fun,x(n)+0.5*h,y(n)+0.5*h

22、*K2); K4=feval(fun,x(n)+h,y(n)+h*K3); y(n+1)=y(n)+h/6*(K1+2*K2+2*K3+K4);endy表1精確解改進Euler公式古典Runge-Kutta公式0.11.0954451150103321.1000000000000001.0954455316930940.21.183215956619923 1.1918181818181821.1832167455059930.31.2649110640673521.2774378337147221.2649122283403920.41.3416407864998741.3582125995

23、602891.3416423537503710.51.4142135623730951.4351329186577961.4142155778900850.61.4832396974191331.5089662535663321.4832422227719920.71.5491933384829671.5803382376552171.5491964523021430.81.6124515496597101.6497834310477111.6124553496589860.91.6733200530681511.7177793478600871.6733246590162561.01.732

24、0508075688771.7847708324979821.7320563651655652.分別用向前Euler公式、向后Euler公式、中心差分公式和梯形公式解常微分方程初值問題取步長。解：M文件： x=0:0.1:1; y=-exp(x);1.應(yīng)用向前Euler法，運用matlab軟件，程序為：h=0.1;x=0:h:1;y=-exp(x);y1=zeros(size(x);n=length(x);y1(1)=1;for i=1:n-1; y1(i+1)=(1+(h/2)*y1(i);enderror=abs(y-y1)2. 應(yīng)用向后Euler法，運用matlab軟件，程序為：h=0.

25、1;x=0:h:1;y=-exp(x);y2=zeros(size(x);n=length(x);y2(1)=1;for i=1:n-1; y2(i+1)=(1/(1-(h/2)*y2(i);enderror=abs(y-y2)3. 應(yīng)用中心差分法，運用matlab軟件，程序為：h=0.1;x=0:h:1;y=-exp(x);y3=zeros(size(x);n=length(x);y3(1)=1;y3(2)=-exp(h);for i=1:n-1; y3(i+2)=y3(i)+h*y3(i+1);enderror=abs(y-y3)4. 應(yīng)用梯形法，運用matlab軟件，程序為：h=0.1;

26、x=0:h:1;y=-exp(x);y4=zeros(size(x);n=length(x);y4(1)=1;for i=1:n-1; y4(i+2)=(8/(4-h)-1)*y4(i);enderror=abs(y-y4).3.分別用一級隱式中點公式、二級隱式中點公式解常微分方程初值問題取步長，計算結(jié)果保存于下表1中，精確解為。解：fun函數(shù)： function f=fun(x,y) f=y-2*x/y;1. 應(yīng)用一級隱式中點公式法，運用matlab軟件，程序為：format longx=0:0.1:1.0;y=sqrt(1+2*x);h=0.1;x=0:h:1.0;y1=zeros(siz

27、e(x);y1=sqrt(1+2*x);y1(1)=1;N=length(x);for n=1:N-1 error=1; K1=feval(fun,x(n),y1(n); while(error>=1.0e-010) K11=feval(fun,x(n)+0.5*h,y1(n)+0.5*h*K1); error=abs(K11-K1); K1=K11; end y1(n+1)=y1(n)+h*K1;ende=abs(y-y1)2. 應(yīng)用二級隱式中點公式法，運用matlab軟件，程序為：format longx=0:0.1:1.0;y=sqrt(1+2*x);h=0.1;x=0:h:1.0

28、;y2=zeros(size(x);y2=sqrt(1+2*x);y2(1)=1;N=length(x);for n=1:N-1 error=1; K1=feval(fun,x(n),y2(n); K2=feval(fun,x(n),y2(n); while(error>=1.0e-010)K11=feval(fun,x(n)+(3-sqrt(3)/6*h,y2(n)+1/4*h*K1+(3-2*sqrt(3)/12*h*K2);K22=feval(fun,x(n)+(3+sqrt(3)/6*h,y2(n)+1/4*h*K2+(3+2*sqrt(3)/12*h*K1); error=ma

29、x(abs(K11-K1),abs(K22-K2); K1=K11; K2=K22; end y2(n+1)=y2(n)+0.5*h*(K1+K2);ende=abs(y-y2)4.分別用Adans顯式公式和Adans隱式公式解常微分方程初值問題取步長。解：fun函數(shù)： function f=fun(x,y); f=-y+x.2;1.應(yīng)用Adans顯式法，運用matlab軟件，程序為：format longh=0.1;x=0:h:1.0;y=zeros(size(x);y1=x.2-2*x+4-3*exp(-x);y(1)=1;N=length(x);for n=1:N-1 K1=feval(

30、fun,x(n),y(n); K2=feval(fun,x(n)+0.5*h,y(n)+0.5*h*K1); K3=feval(fun,x(n)+0.5*h,y(n)+0.5*h*K2); K4=feval(fun,x(n)+h,y(n)+h*K3); y(n+1)=y(n)+h/6*(K1+2*K2+2*K3+K4);endfor n=5:N-1 M1=feval(fun,x(n-1),y(n-1); M2=feval(fun,x(n-2),y(n-2); M3=feval(fun,x(n-3),y(n-3); M4=feval(fun,x(n-4),y(n-4); y(n)=y(n-1)+

31、h/24*(55*M1-59*M2+37*M3-9*M4);endy=abs(y-y1)2.應(yīng)用Adans隱式法，運用matlab軟件，程序為：format longh=0.1;x=0:h:1.0;y=zeros(size(x);y1=x.2-2*x+4-3*exp(-x);y(1)=1;N=length(x);for n=1:N-1 K1=feval(fun,x(n),y(n); K2=feval(fun,x(n)+0.5*h,y(n)+0.5*h*K1); K3=feval(fun,x(n)+0.5*h,y(n)+0.5*h*K2); K4=feval(fun,x(n)+h,y(n)+h*

32、K3); y(n+1)=y(n)+h/6*(K1+2*K2+2*K3+K4);endfor n=3:N-1 error=1 M2=feval(fun,x(n),y(n); M3=feval(fun,x(n-1),y(n-1); M4=feval(fun,x(n-2),y(n-2); temp1=y(n); While(error>=1.0e-013) M1=feval(fun,x(n+1),temp1); tempt2=y(n)+h/24*(9*M1+19*M2-5*M3+M4); error=abs(temp1-temp2); temp1=temp2;endy(n+1)=temp2;e

33、ndy=abs(y-y1).參考文獻1 戴嘉尊，邱建賢，微分方程數(shù)值解.南京：東南大學(xué)出版社，2002.2 馮康. 數(shù)值計算方法.北京:國防工業(yè)出版社，1978.3 周義倉，常微分方程及其應(yīng)用.北京，科學(xué)出版社.2009.4 李瑞遐，何志慶.微分方程數(shù)值解法.上海：華東理工大學(xué)出版社.2005.5 徐樹方，高平，張平文.數(shù)值線性代數(shù).北京：北京大學(xué)出版社.2000.6 余徳浩，湯華中.微分方程數(shù)值解法.北京：高等教育出版社.2003.7 曾金平.微分計算方法.長沙：湖南大學(xué)出版社.2006.8 譚浩強，C程序設(shè)計.北京：清華大學(xué)出版社.2009.9 陸金甫，關(guān)治偏微分方程數(shù)值解法.北京：清華大

34、學(xué)出版社.1987.10 王高雄，周之銘，朱思銘，王壽松.常微分方程.北京：高等教育出版社. 1987.11 胡健偉，湯懷民.微分方程數(shù)值方法.北京:科學(xué)出版社.2000.12 李信真，車剛明，歐陽潔，封建湖.計算方法.西安:西北工業(yè)大學(xué)出版社. 2000.13 瑪康.數(shù)值計算方法.杭州:浙江大學(xué)出版社.2003.14 李榮華，瑪果忱.微分方程數(shù)值方法(第三版).北京:高等教育出版社.1996.15 劉會燈，朱飛.MATLAB編程基礎(chǔ)與典型應(yīng)用.北京:人民郵電出版社.2006.16 蘇金明，阮沈勇.MATLAB實用教程（第二版）.北京:電子工業(yè)出版社.2008.17 劉衛(wèi)國.MATLAB程序設(shè)計教程.北京:中國水利水電出版社.2005.18 馬知恩，周義倉.常微分方程定性與穩(wěn)定性方法.北京:科學(xué)出版社.2007.19 C.W.吉爾.常微分方程初值問題的數(shù)值解法,費景高譯.北京:科學(xué)出版社. 1978.20 樓順天，于衛(wèi)，閏華梁.MATLAB程序設(shè)計語言西安:西安電