第五讲Matlab求解微分方程

第五讲Matlab求解微分方程教学目的：学会用MATLAB求简单微分方程的解析解、数值解，加深对微分方程概念和应用的理解；针对一些具体的问题，如追击问题，掌握利用软件求解微分方程的过程；了解微分方程模型解决问题思维方法及技巧；体会微分方程建摸的艺术性．教学重点：利用机理分析建模微分方程模型，掌握追击问题的建模方法，掌握利用MATLAB求解数值解．教学难点：利用机理分析建模微分方程模型，通过举例，结合图形以及与恰当的假设突破教学难点．1微分方程相关函数（命令）及简介函数名函数功能Dy表示y关于自变量的一阶导数D2y表示y关于自变量的二阶导数dsolve('equ1','equ2',…)求微分方程的解析解，equ1、equ2、…为方程（或条件）simplify(s)对表达式s使用maple的化简规则进行化简[r,how]=simple(s)simple命令就是对表达式s用各种规则进行化简，然后用r返回最简形式，how返回形成这种形式所用的规则．[T,Y]=solver(odefun,tspan,y0)求微分方程的数值解，其中的solver为命令ode45、ode23、ode113、ode15s、ode23s、ode23t、ode23tb之一，odefun是显式常微分方程：00)(),(ytyytfdtdy，在积分区间tspan=],[0ftt上，从0t到ft，用初始条件0y求解，要获得问题在其他指定时间点,210,,ttt上的解，则令tspan=],,,[,210ftttt（要求是单调的）．ezplot(x,y,[tmin,tmax])符号函数的作图命令．x,y为关于参数t的符号函数，[tmin,tmax]为t的取值范围．inline()建立一个内联函数．格式：inline('expr','var1','var2',…)，注意括号里的表达式要加引号．因为没有一种算法可以有效地解决所有的ODE问题，为此，Matlab提供了多种求解器Solver，对于不同的ODE问题，采用不同的Solver．求解器SolverODE类型特点说明ode45非刚性单步算法；4、5阶Runge-Kutta方程；累计截断误差达3)(x大部分场合的首选算法ode23非刚性单步算法；2、3阶Runge-Kutta方程；累计截断误差达3)(x使用于精度较低的情形ode113非刚性多步法；Adams算法；高低精度均可到6310~10计算时间比ode45短ode23t适度刚性采用梯形算法适度刚性情形ode15s刚性多步法；Gear's反向数值微分；精度中等若ode45失效时，可尝试使用ode23s刚性单步法；2阶Rosebrock算法；低精度当精度较低时，计算时间比ode15s短ode23tb刚性梯形算法；低精度当精度较低时，计算时间比ode15s短要特别的是：ode23、ode45是极其常用的用来求解非刚性的标准形式的一阶常微分方程(组)的初值问题的解的Matlab的常用程序，其中：ode23采用龙格-库塔2阶算法，用3阶公式作误差估计来调节步长，具有低等的精度．ode45则采用龙格-库塔4阶算法，用5阶公式作误差估计来调节步长，具有中等的精度．2求解微分方程的一些例子2.1几个可以直接用Matlab求微分方程精确解的例子：例1：求解微分方程22xxexydxdy，并加以验证．求解本问题的Matlab程序为：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))%line4说明：(1)行line1是用命令定义x,y为符号变量．这里可以不写，但为确保正确性，建议写上；(2)行line2是用命令求出的微分方程的解：1/2*exp(-x^2)*x^2+exp(-x^2)*C1(3)行line3使用所求得的解．这里是将解代入原微分方程，结果应该为0，但这里给出：-x^3*exp(-x^2)-2*x*exp(-x^2)*C1+2*x*(1/2*exp(-x^2)*x^2+exp(-x^2)*C1)(4)行line4用simplify()函数对上式进行化简，结果为0，表明)(xyy的确是微分方程的解．例2：求微分方程0'xeyxy在初始条件ey2)1(下的特解，并画出解函数的图形．求解本问题的Matlab程序为：symsxyy=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格式)，即xeeyx，此函数的图形如图1：-6-4-20246-30-20-1001020304050x1/xexp(x)+1/xexp(1)图1y关于x的函数图象2.2用ode23、ode45等求解非刚性的标准形式的一阶常微分方程(组)的初值问题的数值解(近似解)．例３：求解微分方程初值问题1)0(2222yxxydxdy的数值解，求解范围为区间[0,0.5]．fun=inline('-2*y+2*x^2+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.5000y'ans=1.00000.92470.84340.77540.71990.67640.64400.62220.61050.60840.61540.6179图形结果为图2．00.050.10.150.20.250.30.350.40.450.50.60.650.70.750.80.850.90.951图2y关于x的函数图像3常微分在实际中的应用3.1导弹追踪问题1、符号说明，w，乙舰的速率恒为v0；设时刻t乙舰的坐标为((),())XtYt，导弹的坐标为((),())xtyt；当零时刻，((0),(0))(1,0)XY，((0),(0))(0,0)xy．建立微分方程模型．22021(),,01(1)(0)0,'(0)0dykvdydxkxdxxwyy由微分方程模型解得11(1)(1)11(),12(1)2(1)2(1)2(1)kkxxyxkkkkk代入题设的数据1/5k，得到导弹的运行轨迹为4655555(1)(1)81224yxx当1x时245y，即当乙舰航行到点)245,1(处时被导弹击中.被击中时间为:00245vvyt.若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,'*')图３导弹运行轨迹（解析法）图４两种方法对比的导弹运行轨迹2、数值方法求解．设导弹速率恒为w，则得到参数方程为)()()()()()(2222yYyYxXwdtdyxXyYxXwdtdx因乙舰以速度0v沿直线1x运动，设01v，5,1,wXYt，因此导弹运动轨迹的参数方程为：0)0(,0)0()()()1(5)1()()1(52222yxytytxdtdyxytxdtdxMATLAB求解数值解程序如下，结果见图４t0=0,tf=0.21;[t,y]=ode45('eq2',[t0tf],[00]);X=1;Y=0:0.001:0.21;plot(X,Y,'-')plot(y(:,1),y(:,2),'*'),holdonx=0:0.01:1;y=-5*(1-x).^(4/5)/8+5*(1-x).^(6/5)/12+5/24;plot(x,y,'r')很明显数值计算的方法比较简单而且适用．3.2蚂蚁追击问题（1）建立平面直角坐标系．以正多边形的中心为原点,设正多边形的一个顶点为起始点,连接此点和原点作x轴.根据x轴作出相应的y轴;选取足够小的t进行采样．（2）在每一时刻t，计算每只蚂蚁在下一时刻tt时的坐标．不妨设甲追逐对象是乙，在时间t时，甲的坐标为A11(,)xy，乙的坐标为B22(,)xy．甲在tt时在'A点(如图1),其坐标为11(cos,sin)xvtyvt，其中2221212121cos,sin,()()xxyydxxyydd.同理，依次计算下一只蚂蚁在tt时的坐标．通过间隔t进行采样，得到新一轮各个蚂蚁在一个新的正多边形位置坐标．（4）重复2）步，直到d充分小为止．（5）连接每只蚂蚁在各时刻的位置，就得到所求的轨迹．用MALAB求解并作图，函数zhuJi(x,y)在附录一定义，如图６t=[1:8];s=7*exp(t.*2*pi/length(t)*i);x=real(s);y=imag(s);zhuJi(x,y)图6当蚂蚁为7只时的图形习题1.求微分方程0sin2')1(2xxyyx的通解．2.求微分方程xeyyyxsin5'2''的通解．3.求微分方程组'A11,(,)Axy22,(,)Bxyd图1在采样时间内，相连蚂蚁追击00yxdtdyyxdtdx在初始条件0|,1|00ttyx下的特解，并画出解函数()yfx的图形．4.分别用ode23、ode45求上述第3题中的微分方程初值问题的数值解(近似解)，求解区间为[0,2]t．利用画图来比较两种求解器之间的差异．4参考文献:[1]MasteringMatlab6,D.Hanselman,B.Littlefield,清华大学出版,2002[2]赵静等编．数学建模与数学实验（第3版）．北京：高等教育出版社．2008.[3]姜启源编.数学模型（第二版）．北京：高等教育出版社.1993.[4]石勇国.蚂蚁追击问题与等角螺线.宜宾学院学报.2008，(6):23-25.[5]张伟年，杜正东，徐冰．常微分方程．北京：高等教育出版社．2006．5附录附录一:zhuji(x,y)的M文件functionzhuji(x,y)clfv=1;dt=0.05;x(length(x)+1)=x(1);y(length(y)+1)=y(1);plot(x,y,'*')holdforit=1:1000fori=1:length(x)-1d=sqrt((x(i)-x(i+1))^2+(y(i)-y(i+1))^2);x(i)=x(i)+v*dt*(x(i+1)-x(i))/d;y(i)=y(i)+v*dt*(y(i+1)-y(i))/d;endplot(x,y,'.')holdonx(length(x))=x(1);y(length(y))=y(1);end