試驗(yàn)七用matlab求解常微分方程

1、實(shí)驗(yàn)七用matlab求解常微分方程一、實(shí)驗(yàn)?zāi)康模?、熟悉常微分方程的求解方法，了解狀態(tài)方程的概念；2、能熟練使用dsolve函數(shù)求常微分方程(組)的解析解；3、能熟練應(yīng)用ode45ode15s函數(shù)分別求常微分方程的非剛性、剛性的數(shù)值解；4、掌握繪制相圖的方法二、預(yù)備知識(shí)：.微分方程的概念未知的函數(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ù)的最高階

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

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

4、等積分法求解外，大 部分微分方程無(wú)限世界，應(yīng)用中主要依靠數(shù)值解法?？紤]一階常微分方程初值問(wèn)題y(t) f(t, y(t),to t tf y(t0)yo其中 y( yl, y2, ym ) , f ( f1,f2,fm ) ,y0( y10 1y20, ym0 ).所謂數(shù)值解法，就是尋求y在一系列離散節(jié)點(diǎn)t0 t1tn tf上的近似值yk，k 01，n稱hktk 1tk為步長(zhǎng)，通常取為常量h。最簡(jiǎn)單的數(shù)值解法是Euler法。Euler法的思路極其簡(jiǎn)單：在節(jié)點(diǎn)出用差商近似代替導(dǎo)數(shù)y(tky(tk)yD y(tk)h這樣導(dǎo)出計(jì)算公式(稱為Euler格式)Yk1 Yk hf (tk,yk),k0,1

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

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

7、s用來(lái)求解剛性方程組，是用格式同ode45?？梢杂胔elpdsolve, help ode45查閱有關(guān)這些命令的詳細(xì)信息例1求下列微分方程的解析解y ay b y sin(2x) y,y(0) 0,y(0) 1 f f g，g g f，f 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(0)=0,Dy(0)=1,x)simplify(s) %以最簡(jiǎn)形式顯示 s結(jié)果為s =(-1/6*cos(3*

8、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是一一個(gè)結(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ù)值解，并

9、進(jìn)行比較。由clear;s=dsolve(Dy=-y+t+1,y(0)=1,t)simplify(s)可得解析解為y t e t O下面再求其數(shù)值解，先編寫(xiě) M文件fun8.m%MK數(shù) 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; %保留已經(jīng)畫(huà)好的圖形，如果下面再畫(huà)圖，兩個(gè)圖形和并在一起t,y=ode45(fun8,0,1,1);plot(t,y,ro); %畫(huà)數(shù)值解圖形，用紅色小圈畫(huà)xlabel(t),ylabel(y)結(jié)果見(jiàn)圖7.

10、1 TOC o 1-5 h z 1.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解析解與數(shù)值解由圖16.6.1可見(jiàn)，解析解和數(shù)值解吻合得很好。例3求方程ml mg sin , (0)0, (0) 0的數(shù)值解.不妨取l1,g 9.8, (0) 15 .則上面方程可化為9.8sin , (0) 15, (0) 0先看看有沒(méi)有解析解.運(yùn)行MATLAB代碼clear;s=dsolve(D2y=9.8*sin(y),y(0)=15,Dy(0)=0,t)simplify(s)知原方程沒(méi)有解析解.下面求數(shù)值

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向量必須為一列向量運(yùn)行MATLAB代碼clear; close;t,y=ode45(fun9,0,10,15,0);plot(t,y(:,1); %畫(huà)隨時(shí)間變化圖,y(：2) 則表示的值xlabel(t),ylabel(y1)結(jié)果見(jiàn)圖7.20123456t0123456t78910圖7.2數(shù)值解由圖7.2可見(jiàn)， 隨時(shí)間t周期變化。以上這些都是常微分方程的精確解法，也稱為

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

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

14、問(wèn)題在其他指定點(diǎn)to，ti，t2，L上的解，則令tspan tokL，匕(要求 t 單調(diào))，y0初始條件。2例5:求解常微分方程y 2y 2x 2x , o x 0.5, y(0)1的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.

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