實驗七__用matlab求解常微分方程

1、實驗七 用matlab求解常微分方程一、實驗目的： 1、熟悉常微分方程的求解方法，了解狀態(tài)方程的概念； 2、能熟練使用dsolve函數(shù)求常微分方程（組）的解析解； 3、能熟練應用ode45ode15s函數(shù)分別求常微分方程的非剛性、剛性的數(shù)值解； 4、掌握繪制相圖的方法二、預備知識：1微分方程的概念未知的函數(shù)以及它的某些階的導數(shù)連同自變量都由一已知方程聯(lián)系在一起的方程稱為微分方程。如果未知函數(shù)是一元函數(shù)，稱為常微分方程。常微分方程的一般形式為如果未知函數(shù)是多元函數(shù)，成為偏微分方程。聯(lián)系一些未知函數(shù)的一組微分方程組稱為微分方程組。微分方程中出現(xiàn)的未知函數(shù)的導數(shù)的最高階解數(shù)稱為微分方程的階。若方程中

2、未知函數(shù)及其各階導數(shù)都是一次的，稱為線性常微分方程，一般表示為若上式中的系數(shù)均與無關(guān)，稱之為常系數(shù)。2常微分方程的解析解有些微分方程可直接通過積分求解.例如,一解常系數(shù)常微分方程可化為,兩邊積分可得通解為.其中為任意常數(shù).有些常微分方程可用一些技巧,如分離變量法,積分因子法,常數(shù)變異法,降階法等可化為可積分的方程而求得解析解.線性常微分方程的解滿足疊加原理,從而他們的求解可歸結(jié)為求一個特解和相應齊次微分方程的通解.一階變系數(shù)線性微分方程總可用這一思路求得顯式解。高階線性常系數(shù)微分方程可用特征根法求得相應齊次微分方程的基本解，再用常數(shù)變異法求特解。一階常微分方程與高階微分方程可以互化，已給一個階

3、方程設(shè)，可將上式化為一階方程組反過來，在許多情況下，一階微分方程組也可化為高階方程。所以一階微分方程組與高階常微分方程的理論與方法在許多方面是相通的，一階常系數(shù)線性微分方程組也可用特征根法求解。3微分方程的數(shù)值解法除常系數(shù)線性微分方程可用特征根法求解，少數(shù)特殊方程可用初等積分法求解外，大部分微分方程無限世界，應用中主要依靠數(shù)值解法。考慮一階常微分方程初值問題其中所謂數(shù)值解法，就是尋求在一系列離散節(jié)點上的近似值稱為步長，通常取為常量。最簡單的數(shù)值解法是Euler法。Euler法的思路極其簡單：在節(jié)點出用差商近似代替導數(shù)這樣導出計算公式（稱為Euler格式）他能求解各種形式的微分方程。Euler法

4、也稱折線法。Euler方法只有一階精度，改進方法有二階Runge-Kutta法、四階Runge-Kutta法、五階Runge-Kutta-Felhberg法和先行多步法等，這些方法可用于解高階常微分方程（組）初值問題。邊值問題采用不同方法，如差分法、有限元法等。數(shù)值算法的主要缺點是它缺乏物理理解。4解微分方程的MATLAB命令MATLAB中主要用dsolve求符號解析解，ode45,ode23,ode15s求數(shù)值解。  s=dsolve(方程1, 方程2,初始條件1，初始條件2 ,自變量)  用字符串方程表示，自變量缺省值為t。導數(shù)用D表示，2階導數(shù)用D2表示，以

5、此類推。S返回解析解。在方程組情形，s為一個符號結(jié)構(gòu)。tout,yout=ode45(yprime,t0,tf,y0) 采用變步長四階Runge-Kutta法和五階Runge-Kutta-Felhberg法求數(shù)值解，yprime是用以表示f(t,y)的M文件名，t0表示自變量的初始值，tf表示自變量的終值，y0表示初始向量值。輸出向量tout表示節(jié)點(t0,t1, ,tn)T,輸出矩陣yout表示數(shù)值解，每一列對應y的一個分量。若無輸出參數(shù)，則自動作出圖形。ode45是最常用的求解微分方程數(shù)值解的命令，對于剛性方程組不宜采用。ode23與ode45類似，只是精度低一些。ode12s用來求解剛性

6、方程組，是用格式同ode45?？梢杂胔elp dsolve, help ode45查閱有關(guān)這些命令的詳細信息.例1 求下列微分方程的解析解（1）（2）（3）方程（1）求解的MATLAB代碼為：>>clear; >>s=dsolve('Dy=a*y+b')結(jié)果為s =-b/a+exp(a*t)*C1方程（2）求解的MATLAB代碼為：>>clear; >>s=dsolve('D2y=sin(2*x)-y','y(0)=0','Dy(0)=1','x')>>s

7、implify(s)  %以最簡形式顯示s結(jié)果為s =(-1/6*cos(3*x)-1/2*cos(x)*sin(x)+(-1/2*sin(x)+1/6*sin(3*x)*cos(x)+5/3*sin(x)ans =-2/3*sin(x)*cos(x)+5/3*sin(x)方程（3）求解的MATLAB代碼為：>>clear; >>s=dsolve('Df=f+g','Dg=g-f','f(0)=1','g(0)=1')>>simplify(s.f) %s是一個結(jié)構(gòu)>>sim

8、plify(s.g)結(jié)果為ans =exp(t)*cos(t)+exp(t)*sin(t)ans =-exp(t)*sin(t)+exp(t)*cos(t)例 求解微分方程先求解析解，再求數(shù)值解，并進行比較。由>>clear; >>s=dsolve('Dy=-y+t+1','y(0)=1','t')>>simplify(s)可得解析解為。下面再求其數(shù)值解，先編寫M文件fun8.m%M函數(shù)fun8.mfunction f=fun8(t,y)f=-y+t+1;再用命令>>clear; close; t=

9、0:0.1:1;>>y=t+exp(-t); plot(t,y); %化解析解的圖形>>hold on; %保留已經(jīng)畫好的圖形，如果下面再畫圖，兩個圖形和并在一起>>t,y=ode45('fun8',0,1,1); >>plot(t,y,'ro'); %畫數(shù)值解圖形，用紅色小圈畫>>xlabel('t'),ylabel('y')結(jié)果見圖7.1圖16.6.1 解析解與數(shù)值解由圖16.6.1可見，解析解和數(shù)值解吻合得很好。例3 求方程的數(shù)值解.不妨取.則上面方程可化為先看看有

10、沒有解析解.運行MATLAB代碼>>clear; >>s=dsolve('D2y=9.8*sin(y)','y(0)=15','Dy(0)=0','t')>>simplify(s)知原方程沒有解析解.下面求數(shù)值解.令可將原方程化為如下方程組建立M文件fun9.m如下%M文件fun9.mfunction f=fun9(t,y)f=y(2), 9.8*sin(y(1)' %f向量必須為一列向量運行MATLAB代碼>>clear; close; >>t,y=ode45(

11、'fun9',0,10,15,0); >>plot(t,y(:,1); %畫隨時間變化圖,y(:2)則表示的值>>xlabel('t'),ylabel('y1')結(jié)果見圖7.2 圖7.2 數(shù)值解由圖7.2可見，隨時間周期變化。以上這些都是常微分方程的精確解法，也稱為常微分方程的符號解。但是，我們知道，有大量的常微分方程雖然從理論上講，其解是存在的，但我們卻無法求出其解析解，此時，我們需要尋求方程的數(shù)值解，在求常微分方程數(shù)值解方面，MATLAB具有豐富的函數(shù)，我們將其統(tǒng)稱為solver，其一般格式為：T,Y=sol

12、ver(odefun,tspan,y0) 該函數(shù)表示在區(qū)間tspan=t0,tf上，用初始條件y0求解顯式常微分方程solver為命令ode45，ode23，ode113，ode15s，ode23s，ode23t，ode23tb之一，這些命令各有特點。我們列表說明如下：求解器ODE類型特點說明ode45非剛性一步算法，4,5階Runge-Kutta方法累積截斷誤差大部分場合的首選算法ode23非剛性一步算法，2,3階Runge-Kutta方法累積截斷誤差使用于精度較低的情形ode113非剛性多步法，Adams算法，高低精度均可達到計算時間比ode45短ode23t適度剛性采用梯形算法適度剛性情

13、形ode15s剛性多步法，Gears反向數(shù)值積分，精度中等若ode45失效時，可嘗試使用ode23s剛性一步法，2階Rosebrock算法，低精度。當精度較低時，計算時間比ode15s短odefun為顯式常微分方程中的tspan為求解區(qū)間，要獲得問題在其他指定點上的解，則令（要求單調(diào)），y0初始條件。例5：求解常微分方程，的MATLAB程序如下：fun=inline('-2*y+2*x*x+2*x')；x,y=ode23(fun,0,0.5,1)結(jié)果為：x = 0,0.0400,0.0900,0.1400,0.1900,0.2400,0.2900,0.3400,0.3900,0.4400,0.4900,0.5000y = 1.0000,0.9247,0.8434,0.7754,0.7199,0.6764,0.6440,0.6222,0.6105,0.6084,0.6154,0.6179例6：求解常微分方程的解，并畫出解的圖形。分析：這是一個二階非線性方程，用現(xiàn)成的方法均不能求解，但我們可以通過下面的變換，將二階方程化為一階方程組，即可求解。令：，則得到：接著，編寫vdp.m如下：function fy