Matlab求解微分方程

Matlab求解微分方程教学目的：MATLAB摸的艺术性．微分方程相关函数（命令）及简介函数名 函数名 函数功能Dy y关于自变量的一阶导数D2y 表示y关于自变量的二阶导数dsolve('equ1','equ2',…) equ2、…为方程（或条件）simplify(s) 对表达式s使用maple的化简规则进行化简simple命令就是对表达式s用各种规则进行化简，然后[r,how]=simple(s) 用rhow返回形成这种形式所用的规则．求微分方程的数值解，其中的solver为命令ode45、ode23ode15sode23sode23tode23tb[T,Y]=solver(odefun,tspan,y0)odefundtdyf(t,y)，yt)y00在积分区间tspan=[t,t]上，从t 到t ，用初始条件y0 f 0 f 0求解，要获得问题在其他指定时间点t,t,t0 1 2,上的解，tspan=[ttt0 1 2,,t]（要求是单调的．fezplot(x,y,[tmin,tmax]) 符号函数的作图命令．x,y为关于参数t的符号函数，[tmin,tmax]为t的取值范围．inline() 'var1','var2',…)注意括号里的表达式要加引号．因为没有一种算法可以有效地解决所有的ODE问题，为此，MatlabSolverODESolver．求解器求解器SolverODE类型特点说明ode45非刚性非刚性ode23tode15s适度刚性刚性5Runge-Kutta方程；累计截断误差达(x)33Runge-Kutta方程；累计截断误差达(x)3多步法；Adams度均可到103~106采用梯形算法多步法；Gear's反向数值微分；精度中等大部分场合的首选算法ode23非刚性使用于精度较低的情形ode113计算时间比ode45短适度刚性情形若ode45试使用ode23s刚性单步法；2阶Rosebrock算法；当精度较低时，计算时低精度ode23tb刚性梯形算法；低精度间比ode15s短短要特别的是：ode23、ode45是极其常用的用来求解非刚性的标准形式的一阶常微分方程(组)的初值问题的解的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下的特解，并画出解函数的图形．求解本问题的Matlab程序为：symsxyy=dsolve('x*Dy+y-exp(x)=0','y(1)=2*exp(1)','x')ezplot(y)exp(1)(Matlab格式y此函数的图形如图1：

eex，x1/xexp(x)+1/x1/xexp(x)+1/xexp(1)50403020100-10-20-30-6-4-20x246图1y关于x的函数图象用ode23ode45(组)的初值问题的数值解(近似解)．dy

2y2x22x例３：求解微分方程初值问题区间[0,0.5]．

dx 的数值解，求解范围为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关于x的函数图像常微分在实际中的应用导弹追踪问题设位于坐标原点的甲舰向位于x轴上点A（1，0）处的乙舰发射导弹，导弹始终对准乙舰。如果乙舰以最大的速度v0

沿平行于y轴的直线行驶，导弹的速度5v0

，求导弹运行的曲线。当乙舰行驶多远时，导弹将它击中？设导弹在tP(x(t),y(tt乙舰的坐标为X(t),Y(t当零X(0),Y(0))(1,0)(x(0),y(0))(0,0)．PQOPP处的切线，即有tyy'0 1x即 tx)y'y （1）0又根据题意，弧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 当x1时y ，即当乙舰航行到点)处时被导弹击中.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,'*')图３导弹运行轨迹（解析法） 图４两种方法对比的导弹运行轨迹2、数值方法求解．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*')结论：导弹大致在（1,0.2）处击中乙舰。w

dx

dy

w2dt dt由于导弹头始终对准乙舰，故导弹的速度平行于乙舰与导弹头位置的差向量，即dxdtXx2,0 消去t，得到参数方程为(Xx)(Xx)2y)2

dydt

Yyw

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

x1运动，设v0

1w5,Xt，因此导弹运动轨迹的参数方程为：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