用MATLAB解常微分方程

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

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

3、ple(cos(x)A2-sin(x)A2)r = cos(2*x) how = combine4. T,Y = solver( odefun,tspan,y0)求微分方程的數(shù)值解.說明：(1)其中的 solver 為命令 ode45、ode23、ode113> ode15s ode23s、ode23t、 ode23tb 之一.odefun是顯式常微分方程:dy f(t，y)y(to) v。(3)在積分區(qū)間tspan=to，tf 上,從t。到tf ,用初始條件y。求解. 要獲得問題在其他指定時(shí)間點(diǎn)to,ti,t2,上的解，則令tspan=to，tl，t2, ,tf(要求是單調(diào)的).(5)

4、因?yàn)闆]有一種算法可以有效地解決所有的ODE問題，為此，Matlab提供了多種求解器 Solver,對于不同的ODE問題，采用不同的Solver.求解器SolverODE類型特點(diǎn)說明ode45非剛性單步算法；4、5階Runge-Kutta方程；累小雌f誤差達(dá)(x)3大部分場合的首選算 法ode23非剛性單步算法；2、3階Runge-Kutta方程；累小雌f誤差達(dá)(x)3使用于精度較低的情 形ode113非剛性多步法；Adams算法；高低精度均可到10 3 10 6計(jì)算時(shí)間比ode45短ode23t適度W性采用梯形算法適度剛性情形ode15s剛性多步法；Gear's反向數(shù)值微分； 精度中等

5、若ode45失效時(shí)，可 嘗試使用ode23s剛性單步法；2階Rosebrock算法； 低精度當(dāng)精度較低時(shí)，計(jì)算時(shí)間比ode15s短ode23tb剛性梯形算法；低精度當(dāng)精度較低時(shí)，計(jì)算時(shí)間比ode15s短(6)要特別的是：ode23、ode45是極其常用的用來求解非剛性的標(biāo)準(zhǔn)形式的一階常微分方程(組)的初值問題的解的 Matlab的常用程序，其中：ode23采用龍格-庫塔2階算法，用3階公式作誤差估計(jì)來調(diào)節(jié)步長，具有 低等的精度.ode45則采用龍格-庫塔4階算法，用5階公式作誤差估計(jì)來調(diào)節(jié)步長，具 有中等的精度.5. ezplot(x,y,tmin,tmax)：符號函數(shù)的作圖命令.x,y為關(guān)于

6、參數(shù)t的符號函數(shù)，tmin,tmax為t的取值范圍.6. inline()：建立一個(gè)內(nèi)聯(lián)函數(shù).格式：inline('expr','var1','var2',),注意 括號里的表達(dá)式要加引號.例：Q = dblquad(inline('y*sin(x)'), pi, 2*pi, 0, pi)三、實(shí)驗(yàn)內(nèi)容1.幾個(gè)可以直接用Matlab求微分方程精確解的例子:例1:求解微分方程dy 2xy xe x ,并加以驗(yàn)證.dx求解本問題的Matlab程序?yàn)椋簊yms x y%line1y=dsolve('Dy+2*x*y=x*exp(-

7、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*exp(-xA2)*x

8、A2+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程序?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格式)，即 y解函數(shù)的圖形如圖1:1/x exp(x)+1/x exp(1

9、)dx5x例3:求微分方程組 dt dydtx 3y圖1et在初始條件x|t o 1, y |t 0 0下的特解,0并畫出解函數(shù)的圖形.求解本問題的Matlab程序?yàn)?syms x y tx,y=dsolve('Dx+5*x+y=exp(t)','Dy-x-3*y=0','x(0)=1','y(0)=0','t')simple(x);simple(y);ezplot(x,y,0,1.3);axis auto微分方程的特解(式子特別長)以及解函數(shù)的圖形均略.2.用ode2& ode45等求解非剛性的標(biāo)準(zhǔn)形式的

10、一階常微分方程(組)的初值問 題的數(shù)值解(近似解).,“ 、TL、 /+、 口二 dy2y 2x2 2x代m ,例4:求解微分方程初值問題dx2y 2x2x的數(shù)值解，求解范圍為區(qū)y(0) 1間0, 0.5.fun=inline('-2*y+2*xA2+2*x','x',y);x,y=ode23(fun,0,0.5,1);x'；y;plot(x,y,'o-')>> x'ans =0.00000.04000.09000.14000.19000.24000.29000.34000.39000.44000.49000.5000

11、>> y'ans =1.00000.92470.84340.77540.71990.67640.64400.62220.61050.60840.61540.6179圖形結(jié)果為圖2.0.950.90.850.80.750.70.650.600.050.10.150.20.250.30.350.40.450.5圖2例5:求解描述振蕩器的經(jīng)典的Ver der Pol微分方程27.y(1y2)當(dāng) y 0,y(0) 1, y'(0) 0,dtdt分析： 令x1y,x2dx，則 Bx2 ,-dx-(1x12)x2x1.dt dtdt先編寫函數(shù)文件verderpol.m：func

12、tion xprime = verderpol(t,x)global mu;xprime = x(2);mu*(1-x(1)A2)*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.2.52 .1.5 -1,I 1.0.5 .-0 .- 0.5 1- 1 .- 1.5 h,-2 -2.5 |J|0510152025303540圖33.用Euler折線法求解前面講到過，能夠求解的微分方程也是十分有

13、限的.下面介紹用Euler折線 法求微分方程的數(shù)值解（近似解）的方法.Euler折線法求解的基本思想是將微分方程初值問題dx f(x,y), y(x0) v?；梢粋€(gè)代數(shù)方程，即差分方程，主要步驟是用差商y(x h) y(x)替代微商崇,于是:y(Xk h)y(xjf 西，y(Xk),記 xk 1xkh, ykv。y(%)y(xk),從而 yk 1 y(xk h),則有y。xk 1yk 1y(x0), xkh,yk hf (xk, yk).k 0,1,2, ,n 1例6:用Euler折線法求解微分方程初值問題dy dx y(0)2x-2", y的數(shù)值解（步長h取0.4）,求解范圍為區(qū)

14、間0,2. 解：本問題的差分方程為x0,y。1,h 0.4,Xk 1Xkh,k 0,1,2,n 1yk2x hf(Xk,yk)(其中：f(x, y) y ).y相應(yīng)的Matlab程序見附錄1數(shù)據(jù)結(jié)果為：01.00000.40000.80001.20001.60001.40002.12333.11454.45932.00006.3074圖形結(jié)果見圖4:圖4特別說明：本問題可進(jìn)一步利用四階 Runge-Kutta法求解，讀者可將兩個(gè)結(jié) 果在一個(gè)圖中顯示，并和精確值比較，看看哪個(gè)更“精確”？（相應(yīng)的Matlab程 序參見附錄2）.四、自己動手1 .求微分方程(x2 1)y' 2xy sin

15、x 0的通解.2 .求微分方程y'' 2y' 5y exsinx的通解.3 .求微分方程組dx一 x y dt dy出在初始條件xlto 1, y It 0 0下的特解，并畫出解函數(shù)y f（x）的圖形.4 .分別用ode2& ode45求上述第3題中的微分方程初值問題的數(shù)值解（近似解），求解區(qū)間為t 0,2.利用畫圖來比較兩種求解器之間的差異.5 .用Euler折線法求解微分方程初值問題,12x2y y ,yy（0） 1的數(shù)值解（步長h取0.1）,求解范圍為區(qū)間0,2.6 .用四階Runge-Kutta法求解微分方程初值問題x ,y y e cosx,y（0）

16、1的數(shù)值解（步長h取0.1）,求解范圍為區(qū)間0,3.四階 Runge-Kutta法的迭代公式為 （Euler折線法實(shí)為一階Runge-Kutta法）：v。 y(x0),xk 1 xk h,yk1 yk h(L 2L2 2L3 L4) 6L1f(xk,y。k 0,1,2, ,n 1hh、L2f (xk-,yk2L1)hhL3 f (xk -, yk - L2)L4 f X h, yk hL3)相應(yīng)的Matlab程序參見附錄2 .試用該方法求解第5題中的初值問題.7.用。de45方法求上述第6題的常微分方程初值問題的數(shù)值解(近似解), 從而利用畫圖來比較兩者間的差異.五、附錄附錄 1: (fulu1.m) clearf=sym('y+2*x/yA2');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;fo