1、實(shí)驗(yàn)六用matlab求解常微分方程1微分方程的概念未知的函數(shù)以及它的某些階的導(dǎo)數(shù)連同自變量都由一已知方程聯(lián)系在一起的方程稱為 微分方程。如果未知函數(shù)是一元函數(shù)，稱為常微分方程。常微分方程的一般形式為F(t,y,y'y, ,y(n) =0如果未知函數(shù)是多元函數(shù)，成為偏微分方程。聯(lián)系一些未知函數(shù)的一組微分方程組稱為微分方程組。微分方程中出現(xiàn)的未知函數(shù)的導(dǎo)數(shù)的最高階解數(shù)稱為微分方程的階。若方程中未知函數(shù)及其各階導(dǎo)數(shù)都是一次的，稱為線性常微分方程，一般表示為y(n)aQy"an(t)y' an(t)y =b(t)若上式中的系數(shù)ai(t)，i =1,2,，n均與t無關(guān)，稱之為常

2、系數(shù)。2 常微分方程的解析解有些微分方程可直接通過積分求解.例如，一解常系數(shù)常微分方程dydty 1可化為dyy 1=dt,兩邊積分可得通解為t dy=ce -1.其中c為任意常數(shù).有些常微分方程可用一些技巧,如分離變量法，積分因子法，常數(shù)變異法，降階法等可化為可積分的方程而求得解析解線性常微分方程的解滿足疊加原理，從而他們的求解可歸結(jié)為求一個特解和相應(yīng)齊次微 分方程的通解一階變系數(shù)線性微分方程總可用這一思路求得顯式解。高階線性常系數(shù)微分 方程可用特征根法求得相應(yīng)齊次微分方程的基本解，再用常數(shù)變異法求特解。一階常微分方程與高階微分方程可以互化，已給一個n階方程y(n) = f(t，ydy(z)

3、(n d)設(shè)二y，y2二y， ，yn二y，可將上式化為一階方程組M' = y2y2= y3y'n_L = ynM'= f(t，y1，y2,，yn)反過來，在許多情況下，一階微分方程組也可化為高階方程。所以一階微分方程組與高階常微分方程的理論與方法在許多方面是相通的，一階常系數(shù)線性微分方程組也可用特征根法求解。3 微分方程的數(shù)值解法除常系數(shù)線性微分方程可用特征根法求解，少數(shù)特殊方程可用初等積分法求解外，大部分微分方程無限世界，應(yīng)用中主要依靠數(shù)值解法。考慮一階常微分方程初值問題y'(t) = f (t，y(t)，t。£t <tf(to) = yo其中

4、 y = (yi,y2,ym)', f =(fi, f2,fm)',y。=(yio,y2°,ym。)'.所謂數(shù)值解法，就是尋求y(t)在一系列離散節(jié)點(diǎn)to : b tn乞tf上的近似值*,k =o,1，n稱 hk二tk1 -tk為步長，通常取為常量 h。最簡單的數(shù)值解法是 Euler法。Euler法的思路極其簡單：在節(jié)點(diǎn)出用差商近似代替導(dǎo)數(shù)yl(t )y(tk 1) - y(tk)y (tk) h這樣導(dǎo)出計(jì)算公式(稱為 Euler格式)y- =yk hf(tk,yQ,k =0,1,2,他能求解各種形式的微分方程。Euler法也稱折線法。Euler方法只有一階精

5、度，改進(jìn)方法有二階Runge-Kutta法、四階Runge-Kutta法、五階Runge-Kutta-Felhberg法和先行多步法等，這些方法可用于解高階常微分方程(組)初值問 題。邊值問題采用不同方法，如差分法、有限元法等。數(shù)值算法的主要缺點(diǎn)是它缺乏物理理解。4.解微分方程的 MATLAB命令MATLAB中主要用dsolve求符號解析解，ode45,ode23,ode15s求數(shù)值解。s=dsolve(方程1'方程2初始條件1','初始條件2',自變量) 用字符串方程表示，自變量缺省值為t。導(dǎo)數(shù)用D表示，2階導(dǎo)數(shù)用D2表示，以此類推。S返回解析解。在方程組情形

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

7、可以用helpdsolve, help ode45查閱有關(guān)這些命令的詳細(xì)信息例1求下列微分方程的解析解(1)y' = ay b(2) y'sin(2x)-y,y(0) =0,y'(0) =1(3) ff g,gm(0) J,g'(0)方程(1)求解的MATLAB 代碼為：>>clear;>>s=dsolve('Dy=a*y+b')結(jié)果為s =-b/a+exp(a*t)*C1方程(2)求解的MATLAB 代碼為：>>clear;>>s=dsolve('D2y=si n(2*x)-y',

8、'y(0)=0','Dy(0)=1','x')>>simplify(s) %以最簡形式顯示 s結(jié)果為s =(-1/6*cos(3*x)-1/2*cos(x)*si n(x)+(-1/2*si n(x)+1/6*si n(3*x)*cos(x)+5/3*si n(x) ans =-2/3*si n(x)*cos(x)+5/3*si n(x)方程(3)求解的MATLAB 代碼為：>>clear;>>s=dsolve('Df=f+g','Dg=g-f','f(0)=1'

9、,'g(0)=1')>>simplify(s.f) %s是一個結(jié)構(gòu)>>simplify(s.g)結(jié)果為ans =exp(t)*cos(t)+exp(t)*si n(t)ans =-exp(t)*si n( t)+exp(t)*cos(t)例2求解微分方程y'= -y t 1, y(0)=1,先求解析解，再求數(shù)值解，并進(jìn)行比較。由>>clear;>>s=dsolve('Dy=-y+t+1','y(0)=1','t')>>simplify(s)t可得解析解為y = t

10、 e 。下面再求其數(shù)值解，先編寫M文件fun8.m%M®數(shù) fun8.mfun cti on f=fun 8(t,y)f=-y+t+1;再用命令>>clear; close; t=0:0.1:1;>>y=t+exp(-t); plot(t,y); %化解析解的圖形>>hold on; %保留已經(jīng)畫好的圖形，如果下面再畫圖，兩個圖形和并在一起>>t,y=ode45('fu n8',0,1,1);>> plot(t,y,'ro'); %畫數(shù)值解圖形，用紅色小圈畫>>xlabel(

11、9;t'),ylabel('y')結(jié)果見圖7.11.41.35 n1.3 L1.25.y1.2-H1.15-1.1-1.05-1|JLI00.20.40.60.81t圖16.6.1解析解與數(shù)值解由圖16.6.1可見，解析解和數(shù)值解吻合得很好。例3求方程mH" = mgsin 匕二(0)-如"(0) = 0的數(shù)值解.不妨取1 =1，g =9.8(0) =15.則上面方程可化為護(hù)= 9.8sin 二 班0) =15，于(0) =0先看看有沒有解析解.運(yùn)行MATLAB代碼>>clear;>>s=dsolve('D2y=9.8

12、*si n(y)','y(0)=15','Dy(0)=0','t')>>simplify(s)知原方程沒有解析解.下面求數(shù)值解.令 =ly2二F可將原方程化為如下方程組y= y2y2 = 9.8si n(yjy1(0) - 15, y2 (0) = 0建立 M文件fun9.m如下%皿文件fun9.mfun cti on f=fun 9(t,y)f=y(2), 9.8*s in (y(1)' %f向量必須為一列向量運(yùn)行MATLAB代碼>>clear; close;>>t,y=ode45('

13、fu n9',0,10,15,0);>>plot(t,y(:,1); %畫71隨時間變化圖,y(:2)貝U表示“的值>>xlabel('t'),ylabel('y1')結(jié)果見圖7.216.516 - -yIi /i.|15.5 / / 1 / 1115 |iiiiijiii012345678910t圖7.2數(shù)值解由圖7.2可見，二隨時間t周期變化。習(xí)題16-61 .求下列微分方程的解析解F3'=2尹血乂(3) F += 工(。0)(4) yyn-yi2 -1 = 0 j3dx + 2(xa -xyQyfy= 0, 九=1316) / + f + y-cosx,