Matlab求解微分方程

Matlab求解微分方程教學(xué)目的：MATLAB摸的藝術(shù)性．微分方程相關(guān)函數(shù)（命令）及簡介函數(shù)名 函數(shù)名 函數(shù)功能Dy y關(guān)于自變量的一階導(dǎo)數(shù)D2y 表示y關(guān)于自變量的二階導(dǎo)數(shù)dsolve('equ1','equ2',…) equ2、…為方程（或條件）simplify(s) 對表達(dá)式s使用maple的化簡規(guī)則進(jìn)行化簡simple命令就是對表達(dá)式s用各種規(guī)則進(jìn)行化簡，然后[r,how]=simple(s) 用rhow返回形成這種形式所用的規(guī)則．求微分方程的數(shù)值解，其中的solver為命令ode45、ode23ode15sode23sode23tode23tb[T,Y]=solver(odefun,tspan,y0)odefundtdyf(t,y)，yt)y00在積分區(qū)間tspan=[t,t]上，從t 到t ，用初始條件y0 f 0 f 0求解，要獲得問題在其他指定時間點(diǎn)t,t,t0 1 2,上的解，tspan=[ttt0 1 2,,t]（要求是單調(diào)的．fezplot(x,y,[tmin,tmax]) 符號函數(shù)的作圖命令．x,y為關(guān)于參數(shù)t的符號函數(shù)，[tmin,tmax]為t的取值范圍．inline() 'var1','var2',…)注意括號里的表達(dá)式要加引號．因為沒有一種算法可以有效地解決所有的ODE問題，為此，MatlabSolverODESolver．求解器求解器SolverODE類型特點(diǎn)說明ode45非剛性非剛性ode23tode15s適度剛性剛性5Runge-Kutta方程；累計截斷誤差達(dá)(x)33Runge-Kutta方程；累計截斷誤差達(dá)(x)3多步法；Adams度均可到103~106采用梯形算法多步法；Gear's反向數(shù)值微分；精度中等大部分場合的首選算法ode23非剛性使用于精度較低的情形ode113計算時間比ode45短適度剛性情形若ode45試使用ode23s剛性單步法；2階Rosebrock算法；當(dāng)精度較低時，計算時低精度ode23tb剛性梯形算法；低精度間比ode15s短短要特別的是：ode23、ode45是極其常用的用來求解非剛性的標(biāo)準(zhǔn)形式的一階常微分方程(組)的初值問題的解的Matlab的常用程序，其中：ode23采用龍格-2階算法，用3低等的精度．ode45則采用龍格-4階算法，用5有中等的精度．求解微分方程的一些例子Matlab求微分方程精確解的例子：例1：求解微分方程dy2xyxex2，并加以驗證．dxMatlab程序為：symsxy %line1y=dsolve('Dy+2*x*y=x*exp(-x^2)','x') %line2diff(y,x)+2*x*y-x*exp(-x^2) %line3simplify(diff(y,x)+2*x*y-x*exp(-x^2)) %line4line1x,y正確性，建議寫上；line2是用命令求出的微分方程的解：1/2*exp(-x^2)*x^2+exp(-x^2)*C1line3但這里給出：-x^3*exp(-x^2)-2*x*exp(-x^2)*C1+2*x*(1/2*exp(-x^2)*x^2+exp(-x^2)*C1)line4用simplify()0yy(x)的確是微分方程的解．xy'yex02e下的特解，并畫出解函數(shù)的圖形．求解本問題的Matlab程序為：symsxyy=dsolve('x*Dy+y-exp(x)=0','y(1)=2*exp(1)','x')ezplot(y)exp(1)(Matlab格式y(tǒng)此函數(shù)的圖形如圖1：

eex，x1/xexp(x)+1/x1/xexp(x)+1/xexp(1)50403020100-10-20-30-6-4-20x246圖1y關(guān)于x的函數(shù)圖象用ode23ode45(組)的初值問題的數(shù)值解(近似解)．dy

2y2x22x例３：求解微分方程初值問題區(qū)間[0,0.5]．

dx 的數(shù)值解，求解范圍為y()1fun=inline('-2*y+2*x^2+2*x','x','y');[x,y]=ode23(fun,[0,0.5],1); x; plot(x,y,'o-')>>x'ans=0.00000.04000.09000.14000.19000.24000.29000.34000.39000.44000.49000.5000>>y'ans=1.00000.92470.84340.77540.71990.67640.64400.62220.61050.60840.61540.617910.950.90.8510.950.90.850.80.750.70.650.60 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5圖2y關(guān)于x的函數(shù)圖像常微分在實際中的應(yīng)用導(dǎo)彈追蹤問題設(shè)位于坐標(biāo)原點(diǎn)的甲艦向位于x軸上點(diǎn)A（1，0）處的乙艦發(fā)射導(dǎo)彈，導(dǎo)彈始終對準(zhǔn)乙艦。如果乙艦以最大的速度v0

沿平行于y軸的直線行駛，導(dǎo)彈的速度5v0

，求導(dǎo)彈運(yùn)行的曲線。當(dāng)乙艦行駛多遠(yuǎn)時，導(dǎo)彈將它擊中？設(shè)導(dǎo)彈在tP(x(t),y(tt乙艦的坐標(biāo)為X(t),Y(t當(dāng)零X(0),Y(0))(1,0)(x(0),y(0))(0,0)．PQOPP處的切線，即有tyy'0 1x即 tx)y'y （1）0又根據(jù)題意，弧OP的長度是AQ 的5倍，即x0由1（2）消去t整理得

1ydx

t （2）0(1x)y"1 1y'2,5 ,0x15 y(0)0,y'(0)0由微分方程模型解得5 4 5 6 5y(1x)5 (1x)5 8 12 245 當(dāng)x1時y ，即當(dāng)乙艦航行到點(diǎn))處時被導(dǎo)彈擊中.5 24 24為:ty

5 .

=1,則在t=0.21處被擊中.v 24v 00 0利用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)行軌跡（解析法） 圖４兩種方法對比的導(dǎo)彈運(yùn)行軌跡2、數(shù)值方法求解．y1

y,y2

y'，將方程化為方程組，即1x)y''建立m-文件eq1.m

151151y2'

y'y151151y21

/x)functiondy=eq1(x,y)dy=zeros(2,1);dy(1)=y(2);dy(2)=1/5*sqrt(1+y(1)^2)/(1-x);建立主程序ff11.mx0=0,xf=0.9999[x,y]=ode15s('eq1',[x0,xf],[0,0]);plot(x,y(:,1),'b.');holdony=0:0.01:2;plot(1,y,'b*')結(jié)論：導(dǎo)彈大致在（1,0.2）處擊中乙艦。w

dx

dy

w2dt dt由于導(dǎo)彈頭始終對準(zhǔn)乙艦，故導(dǎo)彈的速度平行于乙艦與導(dǎo)彈頭位置的差向量，即dxdtXx2,0 消去t，得到參數(shù)方程為(Xx)(Xx)2y)2

dydt

Yyw

(Xx)dy w(Xx)2(Xx)2y)2dt因乙艦以速度v0

x1運(yùn)動，設(shè)v0

1w5,Xt，因此導(dǎo)彈運(yùn)動軌跡的參數(shù)方程為：x)2x)2(ty)2

5 x)dy

5 (ty)dt x)2(ty)2x(0)0,y(0)0MATLABm-eq2.mFunctiondy=eq2(t,y)dy=zeros(2,1);dy(1)=5*(1-y(1))/sqrt((1-y(1)^2)+(t-y(2))^2)dy(2)=5*(t-y(2))/sqrt((1-y(1)^2)+(t-y(2))^2)建立主程序ff2