(完整版)實驗七用matlab求解常微分方程

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

2、高階解數稱為微分方程的階。若方程中未知函數及其各階導數都是一次的，稱為線性常微分方程，一般表示為y ai(t)y(n1) an i(t)y an(t)y b(t)若上式中的系數ai (t)，i1，2，n均與t無關，稱之為常系數。2 .常微分方程的解析解dy .y 1有些微分方程可直接通過積分求解.例如，一解常系數常彳分方程 出 可化為dydtty 1，兩邊積分可得通解為 y ce 1.其中c為任意常數.有些常微分方程可用一些技巧,如分離變量法，積分因子法，常數變異法，降階法等可化為可積分的方程而求得解析解.線性常微分方程的解滿足疊加原理，從而他們的求解可歸結為求一個特解和相應齊次微 分方程的通

3、解.一階變系數線性微分方程總可用這一思路求得顯式解。高階線性常系數微分 方程可用特征根法求得相應齊次微分方程的基本解，再用常數變異法求特解。一階常微分方程與高階微分方程可以互化，已給一個 n階方程(n)(n 1)、y f(t, y ,y , , y )(n 1)設y1 y，y2 y，yny ,可將上式化為一階方程組yy2y2, y3 yn 1 Vnyn f (t,y1, y2, yn)反過來，在許多情況下，一階微分方程組也可化為高階方程。所以一階微分方程組與高 階常微分方程的理論與方法在許多方面是相通的，一階常系數線性微分方程組也可用特征根法求解。3 .微分方程的數值解法除常系數線性微分方程可

4、用特征根法求解，少數特殊方程可用初等積分法求解外，大 部分微分方程無限世界，應用中主要依靠數值解法?？紤]一階常微分方程初值問題y(t)f(t,y(t),t0 t tfy(to) v。其中 y( yl, y2,ym ) , f ( f1,f2, fm ) ,y0(y10,y20,Vm。).所謂數值解法，就是尋求y在一系列離散節(jié)點tot1 tn tf上的近似值yk，k。，1，n稱hktk 1tk為步長，通常取為常量h。最簡單的數值解法是Euler法。Euler法的思路極其簡單：在節(jié)點出用差商近似代替導數y(tk)y(tki) y(tk)h這樣導出計算公式(稱為Euler格式)ykiyk hf(tk

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

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

7、類似，只是精度低一些。ode12s用來求解剛性方程組，是用格式同ode45?？梢杂胔elpdsolve, help ode45查閱有關這些命令的詳細信息例1求下列微分方程的解析解(1) y ay b y sin(2x) y,y(0) 0,y(0) 1 f f g,g g f,f(0) 1,g(0) 1方程(1)求解的MATLAB 代碼為：clear;s=dsolve(Dy=a*y+b)結果為5 =-b/a+exp(a*t)*C1方程(2)求解的MATLAB 代碼為：clear;s=dsolve(D2y=sin(2*x)-y,y(0)=0,Dy(0)=1,x)simplify(s) %以最簡形式

8、顯示 s結果為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是一一個結構simplify(s.g)結果為ans =exp(t)*cos(t)+exp(t)*sin(t)ans =-exp(t)*sin(t)+exp(t)*cos(t)例2求解微分方程y y t 1

9、, y(0) 1,先求解析解，再求數值解，并進行比較。由clear;s=dsolve(Dy=-y+t+1,y(0)=1,t)simplify(s)可得解析解為y t e t O下面再求其數值解，先編寫 M文件fun8.m%MK數 fun8.mfunction f=fun8(t,y)f=-y+t+1;再用命令clear; close; t=0:0.1:1;y=t+exp(-t); plot(t,y); %化解析解的圖形hold on; %保留已經畫好的圖形，如果下面再畫圖，兩個圖形和并在一起t,y=ode45(fun8,0,1,1);plot(t,y,ro); %畫數值解圖形，用紅色小圈畫xla

10、bel(t),ylabel(y)結果見圖7.11.4_tIIk,1.35.1.3.1.25.y 1.2.r1.15 一/-1.1-.1.05-.1 1c00.20.40.60.81t圖16.6.1解析解與數值解由圖16.6.1可見，解析解和數值解吻合得很好。例3求方程ml mg sin , (0)0, (0) 0的數值解.不妨取l1,g 9.8, (0) 15 .則上面方程可化為9.8sin , (0) 15, (0) 0先看看有沒有解析解.運行MATLAB代碼clear;s=dsolve(D2y=9.8*sin(y),y(0)=15,Dy(0)=0,t)simplify(s)知原方程沒有解析

11、解.下面求數值解.令y1,y2可將原方程化為如下方程組yV2V29.8sin(y1)y(0)15(0)0建立 M文彳fun9.m如下%MC件 fun9.mfunction f=fun9(t,y)f=y(2), 9.8*sin(y(1); %f向量必須為一列向量運行MATLAB代碼clear; close;t,y=ode45(fun9,0,10,15,0);plot(t,y(:,1); %畫隨時間變化圖,y(：2) 則表示的值xlabel(t),ylabel(y1)結果見圖7.20123456t78910圖7.2數值解由圖7.2可見， 隨時間t周期變化。以上這些都是常微分方程的精確解法，也稱為常

12、微分方程的符號解。但是，我們知道，有大量的常微分方程雖然從理論上講，其解是存在的，但我們卻無法求出其解析解，此時，我們需要尋求方程的數值解，在求常微分方程數值解方面， MATLAB具有豐富的函數，我們將其統(tǒng)稱為 solver,其一般格式為：T,Y=solver（odefun,tspan,y0）該函數表示在區(qū)間tspan=t0,tf上，用初始條件y0求解顯式常微分方程y f(t,y)solver 為命令 ode45, ode23, ode113, ode15s ode23s, ode23t, ode23tb 之,這些命令各有特點。我們列表說明如下：求解 器ODE類型特點說明ode45非剛性算法，

13、4,5 階 Runge-Kutta方法累積截斷誤差（x）大部分場合的首選算 法ode23非剛性算法，2,3 階 Runge-Kutta方法累積截斷誤差（x）使用于精度較低的情 形ode113非剛性多步法，Adams算法，高低精度均可達到10 3 10 6計算時間比ode45短ode23t適度剛性采用梯形算法適度剛性情形ode15s剛性多步法，Gear反向 數值積分，精度中等若ode45失效時， 可嘗試使用ode23s剛性2 階 Rosebrock算法，低精度。當精度較低時，計算時間比ode15s短odefun為顯式常微分方程 寸f&y)中的f (t, y)tspan為求解區(qū)間，要獲得問題在其他

14、指定點 3t1,t2,L上的解，則令tspan t0，t1，t2，L，tf(要求 t 單調)，y0初始條件。2例5:求解常微分方程y2y 2x 2x , 0 x 0.5, y 1的MATLAB程序如下：fun=inline(-2*y+2*x*x+2*x) ； x,y=ode23(fun,0,0.5,1)結果為：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.61

15、05,0.6084,0.6154,0.6179例6:求解常微分方程d2y dt2(1 y* y 0,y(0) 1,y(0)0,、，一 的解，并回出解的圖形。分析：這是一個二階非線性方程，用現(xiàn)成的方法均不能求解，但我們可以通 過下面的變換，將二階方程化為一階方程組，即可求解。dyx2-7一.令：x1y, dt ,7,則得到:dx1x2, Xi(0)1dt2 7(1 x12)x2 x1, x2(0)0dt接著，編寫vdp.m如下：function fy=vdp(t,x)fy=x(2);7*(1-x(1)A2)*x(2)-x(1);再編寫m文件sy12_6.m如下：y0=1;0t,x=ode45(vdp,0,40,y0);y=x(:,1);dy=x(:,2);7 / 7plot(t,y,t,dy)三、實驗內容:1 .求下列微分方程的解析解(1) 了rr +2了，-3尸=6一了一3y = 2/皿乂13) 了“十一尸=血工(口 0) y/-/2 -1 = 0了啾