1、第五講Matlab求解微分方程教學(xué)目的：學(xué)會用MATLABt簡單微分方程的解析解、數(shù)值解，加深對微分 方程概念和應(yīng)用的理解；針對一些具體的問題，如追擊問題，掌握利用軟件求解 微分方程的過程；了解微分方程模型解決問題思維方法及技巧；體會微分方程建 摸的藝術(shù)性.教學(xué)重點：利用機(jī)理分析建模微分方程模型， 掌握追擊問題的建模方法，掌 握利用MATLA球:解數(shù)值解.教學(xué)難點：利用機(jī)理分析建模微分方程模型， 通過舉例，結(jié)合圖形以及與恰 當(dāng)?shù)募僭O(shè)突破教學(xué)難點.1微分方程相關(guān)函數(shù)（命令）及簡介函數(shù)名函數(shù)功能Dy表示y關(guān)于自變量的一階導(dǎo)數(shù)D2y表示y關(guān)于自變量的二階導(dǎo)數(shù)dsolve('equ1'

2、,'equ2',)求微分方程的解析解，equl、equ2、為方程（或條件）simplify(s)對表式s使用maple的化簡規(guī)則進(jìn)行化簡r,how=simple(s)simple命令就是對表達(dá)式 s用各種規(guī)則進(jìn)行化 簡，然后用r返回最簡形式，how返回形成這種 形式所用的規(guī)則.求微分方程的數(shù)值解，其中的solver 為命令ode45、ode23、ode113、ode15s、ode23s、ode23t、T,Y=ode23tb之一，odefun是顯式常微分方 程：solver(odefun,tspan,y0)dy f(t y)dt y ,在積分區(qū)間tspan=to，tf 上，從to

3、 y(to) V。至hf,用初始條件y。求解，要獲得問題在其他指定時間點to,ti,t2,上的解，則令tspan=to，ti，t2,tf (要求是單調(diào)的).ezplot(x,y,tmin,tmax)符號函數(shù)的作圖命令.x,y為關(guān)于參數(shù)t的符號函數(shù)，tmin,tmax 為t的取值范圍.inline()建立一個內(nèi)聯(lián)函數(shù).格式：inline('expr','varl', 'var2',)，注意括號里的表達(dá)式要加引號.因為沒有一種算法可以有效地解決所有的ODE問題，為此，Matlab提供了多種求解器Solver ,對于不同的ODE問題，采用不同的Solv

4、er .求解器SolverOD淺型特點說明ode45非剛性單步算法；4、5階Runge-Kutta方程；累小雌f誤差達(dá)(x)3大部分場合的首選算法ode23非剛性單步算法；2、3階Runge-Kutta方程；累小雌f誤差達(dá)(x)3使用于精度較低的情形ode113非剛性多步法；Adams#法；高低精度均可到10 3 10 6計算時間比ode45短ode23t適度W性采用梯形算法適度剛性情形ode15s剛性多步法；Gear's反向數(shù)值微分；精度中等若ode45失效時，可嘗試使用ode23s剛性單步法；2階Rosebrock 算法；低精度當(dāng)精度較低時，計算時間比ode15s短ode23tb剛

5、性梯形算法；低精度當(dāng)精度較低時，計算時間比ode15s短要特別的是：ode23、ode45是極其常用的用來求解非剛性的標(biāo)準(zhǔn)形式的一階常微分方程(組)的初值問題的解的Matlab的常用程序，其中：ode23采用龍格-庫塔2階算法，用3階公式作誤差估計來調(diào)節(jié)步長，具有 低等的精度.ode45則采用龍格-庫塔4階算法，用5階公式作誤差估計來調(diào)節(jié)步長，具 有中等的精度.2求解微分方程的一些例子幾個可以直接用Matlab求微分方程精確解的例子：例1:求解微分方程包2xy xe x ,并加以驗證.dx求解本問題的Matlab程序為：syms x y%line1y=dsolve('Dy+2*x*y=

6、x*exp(-xA2)','x')%line2diff(y,x)+2 *x*y-x*exp(-xA2)%line3simplify(diff(y,x)+2*x*y-x*exp(-xA2) %line4說明：(1)行l(wèi)inel是用命令定義x,y為符號變量.這里可以不寫，但為確 保正確性，建議寫上；(2)行l(wèi)ine2是用命令求出的微分方程的解：1/2*exp(-xA2)*xA2+exp(-xA2)*C1(3)行l(wèi)ine3使用所求得的解.這里是將解代入原微分方程，結(jié)果應(yīng)該為0, 但這里給出：-xA3*exp(-xA2)-2*x*exp(-xA2)*C1+2*x*(1/2*ex

7、p(-xA2)*xA2+exp(-xA2)*C1)(4)行l(wèi)ine4 用simplify。函數(shù)對上式進(jìn)行化簡，結(jié)果為 0 ,表明y y(x)的確是微分方程的解.例2:求微分方程xy y ex 0在初始條件y(1) 2e下的特解，并畫出解函 數(shù)的圖形.求解本問題的Matlab程序為：syms x yy=dsolve('x*Dy+y-exp(x)=0','y(1)=2*exp(1)','x')ezplot(y)微分方程的特解為：y=1/x*exp(x)+1/x* exp (1) (Matlab 格式)，即xe ey e,此函數(shù)的圖形如圖1 : x 1

8、/x exp(x)+1/x exp(1) 50 P-40 .-30 L.20 r.10 .-0 -一- 10 上- 20 ,f- 30 L_lIIIIILa-6-4-20246x圖1 y關(guān)于x的函數(shù)圖象用ode23、ode45等求解非剛性的標(biāo)準(zhǔn)形式的一階常微分方程(組)的初值問題的數(shù)值解(近似解).,、口、口= dy 2v 2x22x 乙心上八/+" .".例3:求解微分方程初值問題dx 2y 2x2x的數(shù)值解，求解范圍為y(0) 1區(qū)間0,x'y'fun=inline('-2*y+2*xA2+2*x','x',y);x,y=

9、ode23(fun,0,1);plot(x,y,'o-')>> x'ans =>> y'ans =圖形結(jié)果為圖2 .0.950.90.850.80.750.70.6500.050.10.150.20.250.30.350.40.450.50.6圖2 y關(guān)于x的函數(shù)圖像3常微分在實際中的應(yīng)用導(dǎo)彈追蹤問題1、符號說明，w,乙艦的速率包為V0；設(shè)時刻t乙艦的坐標(biāo)為(X(t),Y)，導(dǎo)彈的坐標(biāo)為(x(t),y(t)；當(dāng)零時刻，(X(0),Y(0) (1,0) , (x(0), y(0) (0,0),建立微分方程模型.d2y dx2J(dy)2,

10、dx(1 x)y(0) 0,y'(0)V0 c,0 w由微分方程模型解得k 1/、 (1 x)y( x)-2(k 1)k 1(1 x)2( k 1)2( k 1)代入題設(shè)的數(shù)據(jù)k1/5,得到導(dǎo)彈的運(yùn)行軌跡為465(1 x)5* x)-5245即當(dāng)乙艦航行到點（1，五）處時被導(dǎo)彈擊中.被擊中時向為：t工Vo5.右v0=1,則在1二處被擊中.24v°利用MALAB乍圖如圖3 .clear, x=0:1; y=-5*(1-x).A(4/5)/8+5*(1-x).A(6/5)/12+5/24; plot(x,y,'*')11圖3導(dǎo)彈運(yùn)行軌跡（解析法）圖4兩種方法對比的

11、導(dǎo)彈運(yùn)行軌跡2、數(shù)值方法求解.設(shè)導(dǎo)彈速率恒為w,則得到參數(shù)方程為dx dt dy dt(X(X x)2 (Y y)222 (Y(X x)2 (Y y)2x)y)因乙艦以速度v0沿直線x 1運(yùn)動，設(shè)v0 1, w5,X 1,Y t,因此導(dǎo)彈運(yùn)動軌跡的參數(shù)方程為:dx dt dy dt x(0)5(1x)2(t y)25.(1 x)2(t y)20,y(0) 0(1(tx)y)MATLA求解數(shù)值解程序如下，結(jié)果見圖4t0=0,tf=;t,y=ode45( 'eq2' ,t0 tf,0 0);X=1;Y=0:;plot(X,Y,'-')plot(y(:,1),y(:,

12、2),'*' ),hold onx=0:1; y=-5*(1-x)A(4/5)/8+5*(1-x)A(6/5)/12+5/24;plot(x,y, 'r')很明顯數(shù)值計算的方法比較簡單而且適用.螞蟻追擊問題(1)建立平面直角坐標(biāo)系.以正多邊形的中心為原點 ，設(shè)正多邊形的一個 頂點為起始點，連接此點和原點作x軸.根據(jù)x軸作出相應(yīng)的y軸；選取足夠小 的t進(jìn)行采樣.(2)在每一時刻t ,計算每只螞蟻在下一時刻t t時的坐標(biāo).不妨設(shè)甲追 逐對象是乙，在時間t時，甲的坐標(biāo)為A(x1,y1),乙的坐標(biāo)為B(x2,y2) .甲在t t 時在A'點(如圖1),其坐標(biāo)為(

13、為 v t cos , y1 v tsin ),其中 cosxx1,sinyy1,d 7(x2x1PG2yf .同理，依次計算下dd ,一只螞蟻在tt時的坐標(biāo).通過間隔t進(jìn)行采樣，得到新一輪各個螞蟻在一個新的正多邊形位置坐標(biāo).(4)重復(fù)2)步，直到d充分小為止.(5)連接每只螞蟻在各時刻的位置，就得到所求的軌跡.A'A,(xi, yi)B,(x2, y2)圖1在采樣時間內(nèi)，相連螞蟻追擊用MALABO解并作圖，函數(shù)zhuJi(x,y)在附錄一定義，如圖6t=1:8;s=7*exp(t.*2*pi/length*i);x=real(s);y=imag(s);zhuJi(x,y)圖6當(dāng)螞蟻為

14、7只時的圖形習(xí)題1 .求微分方程(x2 1) y' 2xy sin x 。的通解.2 .求微分方程y'' 2y' 5y exsinx的通解.3 .求微分方程組dx 出 蟲 出在初始條件x|t0 1, y |t 0 0下的特解，并畫出解函數(shù)y f(x)的圖形.4 .分別用ode23、ode45求上述第3題中的微分方程初值問題的數(shù)值解（ 近似解 ） ，求解區(qū)間為t 0,2 利用畫圖來比較兩種求解器之間的差異4 參考文獻(xiàn) :1 Mastering Matlab 6, D. Hanselman, B. Littlefield, 清華大學(xué)出版,20022 趙靜等編數(shù)學(xué)建模與數(shù)學(xué)實驗（第 3 版） 北京：高等教育出版社 2008.3 姜啟源編 . 數(shù)學(xué)模型（第二版） 北京：高等教育出版社.1993.4 石勇國 . 螞蟻追擊問題與等角螺線. 宜賓學(xué)院學(xué)報. 2008 ， （6）: 23-25.5 張偉年，杜正東，徐冰常微分方程北京：高等教育出版社20065 附錄 附錄一 :zhuji(x,y) 的M儀件function zhuji(x,y)clfv=1;dt=;x(length(x)+1)=x(1);y(length(y)+1)