用MATLAB解常微分方程

1、實(shí)驗(yàn)四 求微分方程的解一、問(wèn)題背景與實(shí)驗(yàn)?zāi)康膶?shí)際應(yīng)用問(wèn)題通過(guò)數(shù)學(xué)建模所歸納而得到的方程，絕大多數(shù)都是微分方程，真正能得到代數(shù)方程的機(jī)會(huì)很少另一方面，能夠求解的微分方程也是十分有限的，特別是高階方程和偏微分方程（組）這就要求我們必須研究微分方程（組）的解法，既要研究微分方程（組）的解析解法（精確解），更要研究微分方程（組）的數(shù)值解法（近似解）對(duì)微分方程（組）的解析解法(精確解)，Matlab 有專門的函數(shù)可以用，本實(shí)驗(yàn)將作一定的介紹本實(shí)驗(yàn)將主要研究微分方程(組)的數(shù)值解法（近似解），重點(diǎn)介紹 Euler 折線法二、相關(guān)函數(shù)（命令）及簡(jiǎn)介1dsolve('equ1','eq

2、u2',)：Matlab 求微分方程的解析解equ1、equ2、為方程（或條件）寫方程（或條件）時(shí)用 Dy 表示y 關(guān)于自變量的一階導(dǎo)數(shù)，用用 D2y 表示 y 關(guān)于自變量的二階導(dǎo)數(shù)，依此類推2simplify(s)：對(duì)表達(dá)式 s 使用 maple 的化簡(jiǎn)規(guī)則進(jìn)行化簡(jiǎn)例如：syms xsimplify(sin(x)2 + cos(x)2)ans=13r,how=simple(s)：由于 Matlab 提供了多種化簡(jiǎn)規(guī)則，simple 命令就是對(duì)表達(dá)式 s 用各種規(guī)則進(jìn)行化簡(jiǎn)，然后用 r 返回最簡(jiǎn)形式，how 返回形成這種形式所用的規(guī)則例如：syms xr,how=simple(cos(

3、x)2-sin(x)2)r = cos(2*x)how = combine4T,Y = solver(odefun,tspan,y0) 求微分方程的數(shù)值解說(shuō)明：(1) 其中的 solver為命令 ode45、ode23、ode113、ode15s、ode23s、ode23t、ode23tb 之一(2) odefun 是顯式常微分方程：(3) 在積分區(qū)間 tspan=上，從到，用初始條件求解(4) 要獲得問(wèn)題在其他指定時(shí)間點(diǎn)上的解，則令 tspan= （要求是單調(diào)的）(5) 因?yàn)闆](méi)有一種算法可以有效地解決所有的 ODE 問(wèn)題，為此，Matlab 提供了多種求解器 Solver，對(duì)于不同的ODE

4、問(wèn)題，采用不同的Solver求解器SolverODE類型特點(diǎn)說(shuō)明ode45非剛性 單步算法；4、5階Runge-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 短ode23t

5、b剛性梯形算法；低精度當(dāng)精度較低時(shí)，計(jì)算時(shí)間比 ode15s 短(6) 要特別的是：ode23、ode45 是極其常用的用來(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)，具有中等的精度5ezplot(x,y,tmin,tmax)：符號(hào)函數(shù)的作圖命令x,y 為關(guān)于參數(shù)t 的符號(hào)函數(shù)，tmin,tmax 為 t 的取值范圍6inline()：建立一個(gè)內(nèi)聯(lián)函數(shù)格式：inline('expr

6、', 'var1', 'var2',) ，注意括號(hào)里的表達(dá)式要加引號(hào)例：Q = dblquad(inline('y*sin(x)'), pi, 2*pi, 0, pi)三、實(shí)驗(yàn)內(nèi)容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-x*exp(-x2) %line3simplify(diff(

7、y,x)+2*x*y-x*exp(-x2) %line4說(shuō)明：(1) 行l(wèi)ine1是用命令定義x,y為符號(hào)變量這里可以不寫，但為確保正確性，建議寫上；(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é)果為 0， 表明的確是微分方程的解例2：求微分方程在初始條件下的特解

8、，并畫(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例3：求微分方程組在初始條件下的特解，并畫(huà)出解函數(shù)的圖形求解本問(wèn)題的 Matlab 程序?yàn)椋簊yms x y tx,y=dsolve('Dx+5*x+y=exp(t)','Dy-x-3*y=0','x(0)

9、=1','y(0)=0','t')simple(x);simple(y);ezplot(x,y,0,1.3);axis auto微分方程的特解(式子特別長(zhǎng))以及解函數(shù)的圖形均略2. 用ode23、ode45等求解非剛性的標(biāo)準(zhǔn)形式的一階常微分方程(組)的初值問(wèn)題的數(shù)值解(近似解)例4：求解微分方程初值問(wèn)題的數(shù)值解，求解范圍為區(qū)間0, 0.5fun=inline('-2*y+2*x2+2*x','x','y');x,y=ode23(fun,0,0.5,1);x'y'plot(x,y,'o

10、-')>> 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.6105 0.6084 0.6154 0.6179圖形結(jié)果為圖 2圖2例 5：求解描述振蕩器的經(jīng)典的 Ver der Pol 微分方程 分析：令則先編寫函數(shù)文件verderpol.m：function xprime = verd

11、erpol(t,x)global mu;xprime = x(2);mu*(1-x(1)2)*x(2)-x(1);再編寫命令文件vdp1.m：global mu;mu = 7;y0=1;0t,x = ode45('verderpol',0,40,y0);x1=x(:,1);x2=x(:,2);plot(t,x1)圖形結(jié)果為圖3圖33. 用 Euler 折線法求解前面講到過(guò)，能夠求解的微分方程也是十分有限的下面介紹用 Euler 折線法求微分方程的數(shù)值解(近似解)的方法Euler 折線法求解的基本思想是將微分方程初值問(wèn)題化成一個(gè)代數(shù)方程，即差分方程，主要步驟是用差商替代微商,于是

12、：記，從而，則有例 6：用 Euler 折線法求解微分方程初值問(wèn)題的數(shù)值解(步長(zhǎng)h取0.4)，求解范圍為區(qū)間0,2解：本問(wèn)題的差分方程為相應(yīng)的Matlab 程序見(jiàn)附錄 1數(shù)據(jù)結(jié)果為： 0 1.0000 0.4000 1.4000 0.8000 2.1233 1.2000 3.1145 1.6000 4.4593 2.0000 6.3074圖形結(jié)果見(jiàn)圖4：圖4特別說(shuō)明：本問(wèn)題可進(jìn)一步利用四階 Runge-Kutta 法求解，讀者可將兩個(gè)結(jié)果在一個(gè)圖中顯示，并和精確值比較，看看哪個(gè)更“精確”？（相應(yīng)的 Matlab 程序參見(jiàn)附錄 2）四、自己動(dòng)手1. 求微分方程的通解2. 求微分方程的通解3. 求

13、微分方程組 在初始條件下的特解，并畫(huà)出解函數(shù)的圖形4. 分別用 ode23、ode45 求上述第 3 題中的微分方程初值問(wèn)題的數(shù)值解(近似解)，求解區(qū)間為利用畫(huà)圖來(lái)比較兩種求解器之間的差異5. 用 Euler 折線法求解微分方程初值問(wèn)題的數(shù)值解(步長(zhǎng)h取0.1)，求解范圍為區(qū)間0,26. 用四階 Runge-Kutta 法求解微分方程初值問(wèn)題的數(shù)值解(步長(zhǎng)h取0.1)，求解范圍為區(qū)間0,3四階 Runge-Kutta 法的迭代公式為(Euler 折線法實(shí)為一階 Runge-Kutta 法)：相應(yīng)的 Matlab 程序參見(jiàn)附錄 2試用該方法求解第5題中的初值問(wèn)題7. 用 ode45 方法求上述第

14、 6 題的常微分方程初值問(wèn)題的數(shù)值解(近似解)，從而利用畫(huà)圖來(lái)比較兩者間的差異五、附錄附錄 1：(fulu1.m)clearf=sym('y+2*x/y2');a=0;b=2;h=0.4;n=(b-a)/h+1;x=0;y=1;szj=x,y;for i=1:n-1y=y+h*subs(f,'x','y',x,y);x=x+h;szj=szj;x,y;endszjplot(szj(:,1),szj(:,2)附錄 2：(fulu2.m)clearf=sym('y-exp(x)*cos(x)');a=0;b=3;h=0.1;n=(b-a)/h+1;x=0;y=1;szj=x,y;for i=1:n-1l1=subs(f,'x','y'