用matlab求解常微分方程

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

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

3、uot;，y(n1)(n 1)設(shè)y1y，y2y， ，yn y，可將上式化為一階方程組%' y2y2 y y'n 1 ynyn' f(t，y1，y2，yn)反過來，在許多情況下，一階微分方程組也可化為高階方程。所以一階微分方程組與高 階常微分方程的理論與方法在許多方面是相通的，一階常系數(shù)線性微分方程組也可用特征根法求解。3 微分方程的數(shù)值解法除常系數(shù)線性微分方程可用特征根法求解，少數(shù)特殊方程可用初等積分法求解外，大局部微分方程無限世界，應用中主要依靠數(shù)值解法。考慮一階常微分方程初值問題y'(t)f(t，y(t)，t° t tfy(to)yo其中 y(yi

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

5、-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表示，以此類推。S返回解析解。在方程組情形，s為一個符號結(jié)構(gòu)。tout,yout=ode45( ypri

6、me ' ,t采用變步長四階Runge-Kutta 法和五階Runge-Kutta-Felhberg法求數(shù)值解，yprime是用以表示f(t,y)的M文 件名，t0表示自變量的初始值，tf表示自變量的終值，y0表示初始向量值。 輸出向量tout表示節(jié)點(t0,t1,tn)T輸出矩陣yout表示數(shù)值解，每一列對 應y的一個分量。假設(shè)無輸出參數(shù)，那么自動作出圖形。ode45是最常用的求解微分方程數(shù)值解的命令，對于剛性方程組不宜采用。ode23與ode45類似，只是精度低一些。ode12s用來求解剛性方程組，是用格式同ode45?？梢杂胔elpdsolve, help ode45查閱有關(guān)這些

7、命令的詳細信息例1求以下微分方程的解析解們 y' ay b2 y'' sin(2x) y,y(0) 0,y'(0) 13 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=dsolve('D2y=sin(2*x)-y','y(

8、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=g-f','f(0)=1','g(0

9、)=1')>>simplify(s.f) %s 是一個結(jié)構(gòu) >>simplify(s.g)結(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ù)值解，并進行比擬。由>>clear;>>s=dsolve('Dy=-y+t+1','y(0)=1','t') >>simplify(s)t 可得解析解為 y t e 。下面再求其數(shù)值

10、解，先編寫 M 文件%M函數(shù)function f=fun8(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('fun8',0,1,1);>>plot(t,y,'ro'); %畫數(shù)值解圖形，用紅色小圈畫>>xlabel('t'),ylabel('y')結(jié)果

11、見圖1.41.351.31.251.151.11.051t圖解析解與數(shù)值解y 1.2由圖可見，解析解和數(shù)值解吻合得很好。 例3求方程ml "mgsin , (0)0, '(0) 0的數(shù)值解不妨取11,g9.8, (0)15 .那么上面方程可化為MATLAB代碼"9.8si n , (0)15, '(0)0運行先看看有沒有解析解>>clear;>>s=dsolve('D2y=9.8*si n(y)','y(0)=15','Dy(0)=0','t')>>simpl

12、ify(s)知原方程沒有解析解下面求數(shù)值解令可將原方程化為如下方程組y2y 9.8si n(y1)y1(0)15, y2(0)0建立M文件如下%M文件fun cti on f=fun 9(t,y)f=y(2), 9.8*s in (y(1)' %f向量必須為一列向量運行MATLAB代碼>>clear; close;>>t,y=ode45('fu n9',0,10,15,0);>>plot(t,y(:,1); %畫隨時間變化圖,y(：2) 那么表示'的值>>xlabel('t'),ylabel('y1')結(jié)果見圖由圖可見,15067圖7.2數(shù)值解隨時間t周期變化。習題16-61 求以下微分方程的解析解(1) T +召亠一3 = 2產(chǎn)鈕 j*+= sinx (o >0)(4) 2 -1 = 0 h血+ 2(b-W妙=0,=13(5) y + / + y=cosx