Matlab求解微分方程

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

2、;equ2',)求微分方程的解析解，equ1、equ2、為方程（或條件）simplify(s)對(duì)表達(dá)式 s 使用 maple 的化簡(jiǎn)規(guī)則進(jìn)行化簡(jiǎn)r,how=simple(s)simple 命令就是對(duì)表達(dá)式 s 用各種規(guī)則進(jìn)行化簡(jiǎn)，然后用 r 返回最簡(jiǎn)形式，how 返回形成這種形式所用的規(guī)則T,Y = solver(odefun,tspan,y0)求微分方程的數(shù)值解，其中的 solver為命令 ode45、ode23、ode113、ode15s、ode23s、ode23t、ode23tb 之一，odefun 是顯式常微分方程：，在積分區(qū)間 tspan=上，從到，用初始條件求解，要獲得問(wèn)題

3、在其他指定時(shí)間點(diǎn)上的解，則令 tspan=（要求是單調(diào)的）ezplot(x,y,tmin,tmax)符號(hào)函數(shù)的作圖命令x,y 為關(guān)于參數(shù)t 的符號(hào)函數(shù)，tmin,tmax 為 t 的取值范圍inline()建立一個(gè)內(nèi)聯(lián)函數(shù)格式：inline('expr', 'var1', 'var2',) ，注意括號(hào)里的表達(dá)式要加引號(hào)因?yàn)闆](méi)有一種算法可以有效地解決所有的 ODE 問(wèn)題，為此，Matlab 提供了多種求解器 Solver，對(duì)于不同的ODE 問(wèn)題，采用不同的Solver求解器SolverODE類型特點(diǎn)說(shuō)明ode45非剛性單步算法；4、5階Runge-

4、Kutta方程；累計(jì)截?cái)嗾`差達(dá)大部分場(chǎng)合的首選算法ode23非剛性單步算法；2、3階Runge-Kutta方程；累計(jì)截?cái)嗾`差達(dá)使用于精度較低的情形ode113非剛性多步法；Adams算法；高低精度均可到計(jì)算時(shí)間比 ode45 短ode23t適度剛性采用梯形算法適度剛性情形ode15s剛性多步法；Gear's反向數(shù)值微分；精度中等若 ode45 失效時(shí)，可嘗試使用ode23s剛性單步法；2階 Rosebrock 算法；低精度當(dāng)精度較低時(shí)，計(jì)算時(shí)間比 ode15s 短ode23tb剛性梯形算法；低精度當(dāng)精度較低時(shí)，計(jì)算時(shí)間比 ode15s 短要特別的是：ode23、ode45 是極其常用的

5、用來(lái)求解非剛性的標(biāo)準(zhǔn)形式的一階常微分方程(組)的初值問(wèn)題的解的 Matlab 的常用程序，其中：ode23 采用龍格-庫(kù)塔2 階算法，用3 階公式作誤差估計(jì)來(lái)調(diào)節(jié)步長(zhǎng)，具有低等的精度ode45 則采用龍格-庫(kù)塔4 階算法，用5 階公式作誤差估計(jì)來(lái)調(diào)節(jié)步長(zhǎng)，具有中等的精度2 求解微分方程的一些例子2.1 幾個(gè)可以直接用 Matlab 求微分方程精確解的例子：例1：求解微分方程，并加以驗(yàn)證求解本問(wèn)題的Matlab 程序?yàn)椋簊yms x y %line1y=dsolve('Dy+2*x*y=x*exp(-x2)','x') %line2diff(y,x)+2 *x*y

6、-x*exp(-x2) %line3simplify(diff(y,x)+2*x*y-x*exp(-x2) %line4說(shuō)明：(1) 行l(wèi)ine1是用命令定義x,y為符號(hào)變量這里可以不寫(xiě)，但為確保正確性，建議寫(xiě)上；(2) 行l(wèi)ine2是用命令求出的微分方程的解：1/2*exp(-x2)*x2+exp(-x2)*C1(3) 行l(wèi)ine3使用所求得的解這里是將解代入原微分方程，結(jié)果應(yīng)該為0，但這里給出：-x3*exp(-x2)-2*x*exp(-x2)*C1+2*x*(1/2*exp(-x2)*x2+exp(-x2)*C1)(4) 行l(wèi)ine4 用 simplify() 函數(shù)對(duì)上式進(jìn)行化簡(jiǎn)，結(jié)果為

7、 0， 表明的確是微分方程的解例2：求微分方程在初始條件下的特解，并畫(huà)出解函數(shù)的圖形求解本問(wèn)題的 Matlab 程序?yàn)椋簊yms 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格式)，即，此函數(shù)的圖形如圖 1：圖1 y關(guān)于x的函數(shù)圖象2.2 用ode23、ode45等求解非剛性的標(biāo)準(zhǔn)形式的一階常微分方程(組)的初值問(wèn)題的數(shù)值解(近似解)例：求解微分方程初值問(wèn)題的數(shù)值解，求解范圍為區(qū)間0, 0.

8、5fun=inline('-2*y+2*x2+2*x','x','y');x,y=ode23(fun,0,0.5,1);x'y'plot(x,y,'o-')>> x'ans =0.0000 0.0400 0.0900 0.1400 0.1900 0.24000.2900 0.3400 0.3900 0.4400 0.4900 0.5000>> y'ans =1.0000 0.9247 0.8434 0.7754 0.7199 0.67640.6440 0.6222 0.610

9、5 0.6084 0.6154 0.6179圖形結(jié)果為圖 2圖2 y關(guān)于x的函數(shù)圖像3 常微分在實(shí)際中的應(yīng)用 3.1 導(dǎo)彈追蹤問(wèn)題1、符號(hào)說(shuō)明，乙艦的速率恒為v0；設(shè)時(shí)刻乙艦的坐標(biāo)為，導(dǎo)彈的坐標(biāo)為；當(dāng)零時(shí)刻，建立微分方程模型由微分方程模型解得代入題設(shè)的數(shù)據(jù)，得到導(dǎo)彈的運(yùn)行軌跡為當(dāng)時(shí)，即當(dāng)乙艦航行到點(diǎn)處時(shí)被導(dǎo)彈擊中. 被擊中時(shí)間為:. 若v0=1, 則在t=0.21處被擊中.利用MALAB作圖如圖clear, x=0:0.01:1; y=-5*(1-x).(4/5)/8+5*(1-x).(6/5)/12+5/24; plot(x,y,'*')圖導(dǎo)彈運(yùn)行軌跡（解析法） 圖兩種方法對(duì)

10、比的導(dǎo)彈運(yùn)行軌跡2、數(shù)值方法求解設(shè)導(dǎo)彈速率恒為，則得到參數(shù)方程為 因乙艦以速度沿直線運(yùn)動(dòng)，設(shè)，因此導(dǎo)彈運(yùn)動(dòng)軌跡的參數(shù)方程為： MATLAB求解數(shù)值解程序如下，結(jié)果見(jiàn)圖t0=0,tf=0.21;t,y=ode45('eq2',t0 tf,0 0);X=1;Y=0:0.001:0.21;plot(X,Y,'-')plot(y(:,1),y(:,2),'*'),hold onx=0:0.01:1; y=-5*(1-x).(4/5)/8+5*(1-x).(6/5)/12+5/24;plot(x,y,'r')很明顯數(shù)值計(jì)算的方法比較簡(jiǎn)單而且

11、適用 3.2 螞蟻?zhàn)窊魡?wèn)題（1）建立平面直角坐標(biāo)系以正多邊形的中心為原點(diǎn), 設(shè)正多邊形的一個(gè)頂點(diǎn)為起始點(diǎn), 連接此點(diǎn)和原點(diǎn)作軸. 根據(jù)軸作出相應(yīng)的軸; 選取足夠小的進(jìn)行采樣（2）在每一時(shí)刻，計(jì)算每只螞蟻在下一時(shí)刻時(shí)的坐標(biāo)不妨設(shè)甲追逐對(duì)象是乙，在時(shí)間時(shí)，甲的坐標(biāo)為，乙的坐標(biāo)為甲在時(shí)在點(diǎn)(如圖1), 其坐標(biāo)為，其中. 同理，依次計(jì)算下一只螞蟻在時(shí)的坐標(biāo)通過(guò)間隔進(jìn)行采樣，得到新一輪各個(gè)螞蟻在一個(gè)新的正多邊形位置坐標(biāo)（4）重復(fù)2）步，直到充分小為止（5）連接每只螞蟻在各時(shí)刻的位置，就得到所求的軌跡圖1 在采樣時(shí)間內(nèi)，相連螞蟻?zhàn)窊粲肕ALAB求解并作圖，函數(shù)zhuJi(x,y)在附錄一定義，如圖t=1:

12、8; s=7*exp(t.*2*pi/length(t)*i); x=real(s); y=imag(s);zhuJi(x,y) 圖6當(dāng)螞蟻為7只時(shí)的圖形習(xí)題1. 求微分方程的通解2. 求微分方程的通解3. 求微分方程組 在初始條件下的特解，并畫(huà)出解函數(shù)的圖形4. 分別用 ode23、ode45 求上述第 3 題中的微分方程初值問(wèn)題的數(shù)值解(近似解)，求解區(qū)間為利用畫(huà)圖來(lái)比較兩種求解器之間的差異4 參考文獻(xiàn):1 Mastering Matlab 6, D. Hanselman, B. Littlefield, 清華大學(xué)出版,20022 趙靜等編數(shù)學(xué)建模與數(shù)學(xué)實(shí)驗(yàn)（第3版）北京：高等教育出版社2008.3 姜啟源編. 數(shù)學(xué)模型（第二版）北京：高等教育出版社.1993.4 石勇國(guó). 螞蟻?zhàn)窊魡?wèn)題與等角螺線. 宜賓學(xué)院學(xué)報(bào). 2008，(6): 23-25.5 張偉年，杜正東，徐冰常微分方程北京：高等教育出版社20065 附錄附錄一:zhuji(x,y)的M文件function zhuji(x,y)clfv=1;dt=0.05;x(length(x)+1)=x(1);y(length(y)+1)=y(1);plot(x,y,'*') holdfor