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

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

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

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

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

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

7、類似，只是精度低一些。ode12s用來求解剛性方程組，是用格式同ode45?？梢杂胔elpdsolve, help ode45查閱有關(guān)這些命令的詳細(xì)信息例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)結(jié)果為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結(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)=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

9、, y(0) 1,先求解析解，再求數(shù)值解，并進(jìn)行比較。由clear;s=dsolve(Dy=-y+t+1,y(0)=1,t)simplify(s)可得解析解為y t e t O下面再求其數(shù)值解，先編寫 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)畫好的圖形，如果下面再畫圖，兩個圖形和并在一起t,y=ode45(fun8,0,1,1);plot(t,y,ro); %畫數(shù)值解圖形，用紅色小圈畫xla

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

11、解.下面求數(shù)值解.令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); %畫隨時間變化圖,y(：2) 則表示的值xlabel(t),ylabel(y1)結(jié)果見圖7.20123456t78910圖7.2數(shù)值解由圖7.2可見， 隨時間t周期變化。以上這些都是常微分方程的精確解法，也稱為常

12、微分方程的符號解。但是，我們知道，有大量的常微分方程雖然從理論上講，其解是存在的，但我們卻無法求出其解析解，此時，我們需要尋求方程的數(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, ode15s ode23s, ode23t, ode23tb 之,這些命令各有特點(diǎn)。我們列表說明如下：求解 器ODE類型特點(diǎn)說明ode45非剛性算法，

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

14、指定點(diǎn) 3t1,t2,L上的解，則令tspan t0，t1，t2，L，tf(要求 t 單調(diào))，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)結(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.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)三、實(shí)驗(yàn)內(nèi)容:1 .求下列微分方程的解析解(1) 了rr +2了，-3尸=6一了一3y = 2/皿乂13) 了“十一尸=血工(口 0) y/-/2 -1 = 0了啾