用matlab求解常微分方程

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

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

3、(n 1)、y f(t, y',y", , y )(n 1)設(shè)y1y,y2y , ,yn y ,可將上式化為一階方程組yy2V2' y3y'n 1 VnYn' f (t,y1, y2, yn)反過來，在許多情況下，一階微分方程組也可化為高階方程。所以一階微分方程組與高 階常微分方程的理論與方法在許多方面是相通的，一階常系數(shù)線性微分方程組也可用特征根法求解。3.微分方程的數(shù)值解法除常系數(shù)線性微分方程可用特征根法求解，少數(shù)特殊方程可用初等積分法求解外，大部分微分方程無限世界，應(yīng)用中主要依靠數(shù)值解法?？紤]一階常微分方程初值問題y'(t)f(t,y(t

4、),t。t tfy(to)yo其中 y(y1,y2,ym )', f (f1,f2, fm)', y0(y10,y20,ym。)'.所謂數(shù)值解法，就是尋求在一系列離散節(jié)點(diǎn)t0t1tn tf上的近似值yk,k0,1，,n稱hk tk 1 tk為步長，通常取為常量h。最簡單的數(shù)值解法是Euler法。Euler法的思路極其簡單：在節(jié)點(diǎn)出用差商近似代替導(dǎo)數(shù)y'(tk)y(tki) y(tk)h這樣導(dǎo)出計(jì)算公式(稱為Euler格式)ykiyk hf(tk,yk),k 0,1,2,他能求解各種形式的微分方程。Euler法也稱折線法。Euler方法只有一階精度，改進(jìn)方法有二階

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

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

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

8、olve('D2y=sin(2*x)-y','y(0)=0','Dy(0)=1','x')>>simplify(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=

9、g-f','f(0)=1','g(0)=1')>>simplify %s 是一個(gè)結(jié)構(gòu)>>simplify結(jié)果為ans =exp(t)*cos(t)+exp(t)*sin(t)ans =-exp(t)*sin(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

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

11、'y')結(jié)果見圖1.41.351.3 .'.1.25ry 1.2-1.15-1.1-1.05-1 L11ei00.20.40.60.81t圖16.6.1解析解與數(shù)值解由圖16.6.1可見，解析解和數(shù)值解吻合得很好。例3求方程ml " mg sin , (0)0, '(0) 0的數(shù)值解.不妨取l 1,g9.8, (0) 15 .則上面方程可化為"9.8sin , (0) 15, '(0) 0先看看有沒有解析解運(yùn)行MATLAB代碼>>clear;>>s=dsolve('D2y=*sin(y)',&#

12、39;y(0)=15','Dy(0)=0','t')>>simplify(s)知原方程沒有解析解.下面求數(shù)值解.令y1,y2'可將原方程化為如下方程組yV2V2'9.8sin(y1)y(0)15,y2(0)0建立M文件如下%M文件function f=fun9(t,y)f=y(2), *sin(y(1)' %f向量必須為一列向量 運(yùn)行MATLAB代碼>>clear; close; >>t,y=ode45('fun9',0,10,15,0);>>plot(t,y(:,1);%畫隨時(shí)間變化圖，y(:2)則表示的值>>xlabel('t'),ylabel('y1')結(jié)果見圖16.51615.534數(shù)值解圖5678910t15 012由圖可見，隨時(shí)間t周期變化。習(xí)題16-61.求下列微分方程的解析解(1) 了" + 2y'-3丁=名了"一3尸=2產(chǎn)皿乂(3)= sin a (q >0) J7ff-yf2 -1 = 0 V必+ 21-呼”沖=0, yx i =13(幻 V + y'+y = ss工，兒一0=。, 了G1/ + y = V+3§瑞 jLL y'L