常微分方程組的MATLAB求解方法

1、一、常微分方程組(ODEs)簡介 11. 簡諧振動12 電路 Vander Pol方程 13. 生物種群的Volterra-Lotka方程 24. 蝴蝶效應Lorenz方程2二、MATLAB數(shù)值求解0DEs的方法 31. 多變量常微分方程組的求解42. 高階常微分方程如何表示？ 43. 相圖和極限環(huán)怎么繪制？ 4個人在學習自動控制原理、現(xiàn)代控制理論、非線性動力學等課程時，經(jīng)常遇到求解常微分方程組的問 題。很多人知道MATLAB是簡便易行的一個工具，但是不會調(diào)用它自帶的ode求解器，往往還在自己 編寫單步歐拉法的程序，不僅求解精度差，而且程序不規(guī)范，還浪費了大量時間。以下我就工程中常 見的一些非

2、線性系統(tǒng)，利用MATLAB自帶的求解器，說明一下如何求解ODE方程組、以及如何繪制 相軌跡和極限環(huán)的問題。供相尖專業(yè)工科大學生參考和借鑒。一、常微分方程組(ODEs)簡介以下列出了一些較為著名的非線性動力學系統(tǒng)的數(shù)學表達式，大都是由常微分方程組表達的。這種形式 在工程中應用非常廣泛，如力學中的非線性振動、航天領域的彈道計算、控制工程中的非線性系統(tǒng)等， 由于自然界的大多數(shù)現(xiàn)象都表現(xiàn)出非線性，因此對于該種動力系統(tǒng)的研究以及微分方程的求解也具有重大 的意義。以下列出一些工程應用中常見的一些由ODE方程組所描述的動力系統(tǒng)。1. 簡諧振動該式是一個2階非線性常微分方程。2. 電路Vander Pol方程

3、2Fig 1VanderPol系統(tǒng)時域響應3. 生物種群的Volterra丄otka方程Fig 2. Volterra丄otka方程時域響應（左）Fig 3非線性動力學方程的極限環(huán)（右）左圖的捕食者獵物隨時間變化的曲線表現(xiàn)出強烈的非線性，而狀態(tài)變量x、y的變化卻呈現(xiàn)出一個規(guī)則的鵝卵石狀。4. 蝴蝶效應Lorenz方程01 0203040 SO 60604020x020-40Fig 4.洛倫茲方程系統(tǒng)3維狀態(tài)的時域響應(左亠s小“亠亠、一Fig 5系統(tǒng)的3維空間相圖：蝴蝶效應(右)世微分方程組的魅力就在于，它可以表現(xiàn)出使人意想不到的規(guī)律性與自然之美。按照數(shù)值積分的方法對3個狀態(tài) 變量xyz進行求

4、解，由于該動力系統(tǒng)較為敏感，因此需要以盡可能小的積分步長來求 解。求解所得的相圖在3維呈現(xiàn)出一個蝴蝶狀。MATLAB數(shù)值求解ODEs的方法對于單個的常微分方程，給定初值和求解區(qū)間，同學們一定會求解。如以下形式的一階非 線性微分方程，典型的方法有：單步Euler法、4階Runge-Kutta積分法等。大多數(shù) 的數(shù)值分析或計算方法教材上均有講述。dyf ( x，y)dxy (xo) yo在MATLAB中，編寫被積函數(shù)f(x,y)的腳本，通過調(diào)用命令ode45(f,yO,xO,xf)可以得到較好的數(shù)值解。然而，難倒大多數(shù)初學者的情況卻是：多變量微分方程，狀態(tài)如何表示，函數(shù)的輸入和輸出該怎么用？單變量

5、二階或高階微分方程，只能求解一階ODE的函數(shù)命令如何用在高階系統(tǒng)上？。要求繪制非線性動力學系統(tǒng)的相軌跡圖，該如何繪制？以下我就這兩個問題進行解釋1-多變量常微分方程組的求解多變量常微分方程組的求解，在MATLAB中仍然可 以按照函數(shù)ode模板進行求解，其形式與一階單變量是一樣的。由于MATLAB的ode 函數(shù)只能接受形參表形式為(t, x)的被積函數(shù)，所以我們要將狀態(tài)變量通過多維矢量x傳 遞進去。如14所述的生物種群2變量微分方程組Valter Lotka，將方程中的狀態(tài)變量X與y分別表示為一個二維向量的2維，從而按照類似方法進行求解。function dx=ValterraLotka(t,x

6、) A=2/3;B=4/3;C=1;D=1; dx=x(1)*(A-B*x(2); dx(2)=x(2)*(C*x(1)-D); dx=dx(:);end2. 高階常微分方程如何表示？對于第2種形式的ODE，般來說我們可以引入一些附加變量而將其化為第1種形式，從而通過第1種形式的求解而得到高階微分的結(jié)果。如Vander Pol方程是一個 二階常微分方程，我們可以將其化為：在MATLAB中，同樣可以編寫多變量的被積函數(shù)VanderPol(t,x)，需要注意的是這里 的狀態(tài)變量x是以上我們降階之后的，即為矢量形式。同樣函數(shù)輸出的也是2變量的矢量形式，即.function dx=VanderPol(

7、t，x)%經(jīng)過測試可得mu影響極限環(huán)的形狀；初值不影響最終形狀mu=0.8;dx(1)=x(2);dx(2)=mu*(1 -x(1)*x(1 )*x(2)-x;dx=dx(:);end3. 相圖和極限環(huán)怎么繪制？對于該二階單變量微分方程，控制工程中有極限環(huán)的概念。 它的定義是：在時間t趨于正無窮時，相軌跡都逐漸靠近某一條封閉的曲線，這種被稱為 穩(wěn)定的極限環(huán)。繪制極限環(huán)或者相圖，首先得求解ODE方程組，得到狀態(tài)的時域響應。 在MATLAB中，極限環(huán)的繪制方法很簡單，因為以上降階的過程中有：所以對微分方程 組的兩個狀態(tài)變量x1和x2是函數(shù)輸出的結(jié)果，對其直接作plot即可得到相平面上狀 態(tài)變量的運動軌跡，即極限環(huán)。t,y=ode45(VanderPol,0,15,-0.5,0); figure(1 ),plot(t,y(:,1 ),t,y(:,2),xlabel(f ),ylabel(刈; figure(2),plot(y(:,1),y(:,2),axis equal, xlabel(x),ylabel(dx/