Matlab求解微分方程_第1页
Matlab求解微分方程_第2页
Matlab求解微分方程_第3页
Matlab求解微分方程_第4页
Matlab求解微分方程_第5页
已阅读5页,还剩2页未读 继续免费阅读

下载本文档

版权说明:本文档由用户提供并上传,收益归属内容提供方,若内容存在侵权,请进行举报或认领

文档简介

1、第五讲 Matlab求解微分方程教学目的:学会用MATLAB求简单微分方程的解析解、数值解,加深对微分方程概念和应用的理解;针对一些具体的问题,如追击问题,掌握利用软件求解微分方程的过程;了解微分方程模型解决问题思维方法及技巧;体会微分方程建摸的艺术性教学重点:利用机理分析建模微分方程模型,掌握追击问题的建模方法,掌握利用MATLAB求解数值解教学难点:利用机理分析建模微分方程模型,通过举例,结合图形以及与恰当的假设突破教学难点1 微分方程相关函数(命令)及简介函数名函数功能Dy表示y 关于自变量的一阶导数D2y表示 y 关于自变量的二阶导数dsolve('equ1','

2、;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 是显式常微分方程:,在积分区间 tspan=上,从到,用初始条件求解,要获得问题

3、在其他指定时间点上的解,则令 tspan=(要求是单调的)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-

4、Kutta方程;累计截断误差达大部分场合的首选算法ode23非刚性单步算法;2、3阶Runge-Kutta方程;累计截断误差达使用于精度较低的情形ode113非刚性多步法;Adams算法;高低精度均可到计算时间比 ode45 短ode23t适度刚性采用梯形算法适度刚性情形ode15s刚性多步法;Gear's反向数值微分;精度中等若 ode45 失效时,可尝试使用ode23s刚性单步法;2阶 Rosebrock 算法;低精度当精度较低时,计算时间比 ode15s 短ode23tb刚性梯形算法;低精度当精度较低时,计算时间比 ode15s 短要特别的是:ode23、ode45 是极其常用的

5、用来求解非刚性的标准形式的一阶常微分方程(组)的初值问题的解的 Matlab 的常用程序,其中:ode23 采用龙格-库塔2 阶算法,用3 阶公式作误差估计来调节步长,具有低等的精度ode45 则采用龙格-库塔4 阶算法,用5 阶公式作误差估计来调节步长,具有中等的精度2 求解微分方程的一些例子2.1 几个可以直接用 Matlab 求微分方程精确解的例子:例1:求解微分方程,并加以验证求解本问题的Matlab 程序为:syms x y %line1y=dsolve('Dy+2*x*y=x*exp(-x2)','x') %line2diff(y,x)+2 *x*y

6、-x*exp(-x2) %line3simplify(diff(y,x)+2*x*y-x*exp(-x2) %line4说明:(1) 行line1是用命令定义x,y为符号变量这里可以不写,但为确保正确性,建议写上;(2) 行line2是用命令求出的微分方程的解:1/2*exp(-x2)*x2+exp(-x2)*C1(3) 行line3使用所求得的解这里是将解代入原微分方程,结果应该为0,但这里给出:-x3*exp(-x2)-2*x*exp(-x2)*C1+2*x*(1/2*exp(-x2)*x2+exp(-x2)*C1)(4) 行line4 用 simplify() 函数对上式进行化简,结果为

7、 0, 表明的确是微分方程的解例2:求微分方程在初始条件下的特解,并画出解函数的图形求解本问题的 Matlab 程序为:syms 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格式),即,此函数的图形如图 1:图1 y关于x的函数图象2.2 用ode23、ode45等求解非刚性的标准形式的一阶常微分方程(组)的初值问题的数值解(近似解)例:求解微分方程初值问题的数值解,求解范围为区间0, 0.

8、5fun=inline('-2*y+2*x2+2*x','x','y');x,y=ode23(fun,0,0.5,1);x'y'plot(x,y,'o-')>> 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.610

9、5 0.6084 0.6154 0.6179图形结果为图 2图2 y关于x的函数图像3 常微分在实际中的应用 3.1 导弹追踪问题1、符号说明,乙舰的速率恒为v0;设时刻乙舰的坐标为,导弹的坐标为;当零时刻,建立微分方程模型由微分方程模型解得代入题设的数据,得到导弹的运行轨迹为当时,即当乙舰航行到点处时被导弹击中. 被击中时间为:. 若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,'*')图导弹运行轨迹(解析法) 图两种方法对

10、比的导弹运行轨迹2、数值方法求解设导弹速率恒为,则得到参数方程为 因乙舰以速度沿直线运动,设,因此导弹运动轨迹的参数方程为: MATLAB求解数值解程序如下,结果见图t0=0,tf=0.21;t,y=ode45('eq2',t0 tf,0 0);X=1;Y=0:0.001:0.21;plot(X,Y,'-')plot(y(:,1),y(:,2),'*'),hold onx=0:0.01:1; y=-5*(1-x).(4/5)/8+5*(1-x).(6/5)/12+5/24;plot(x,y,'r')很明显数值计算的方法比较简单而且

11、适用 3.2 蚂蚁追击问题(1)建立平面直角坐标系以正多边形的中心为原点, 设正多边形的一个顶点为起始点, 连接此点和原点作轴. 根据轴作出相应的轴; 选取足够小的进行采样(2)在每一时刻,计算每只蚂蚁在下一时刻时的坐标不妨设甲追逐对象是乙,在时间时,甲的坐标为,乙的坐标为甲在时在点(如图1), 其坐标为,其中. 同理,依次计算下一只蚂蚁在时的坐标通过间隔进行采样,得到新一轮各个蚂蚁在一个新的正多边形位置坐标(4)重复2)步,直到充分小为止(5)连接每只蚂蚁在各时刻的位置,就得到所求的轨迹图1 在采样时间内,相连蚂蚁追击用MALAB求解并作图,函数zhuJi(x,y)在附录一定义,如图t=1:

12、8; s=7*exp(t.*2*pi/length(t)*i); x=real(s); y=imag(s);zhuJi(x,y) 图6当蚂蚁为7只时的图形习题1. 求微分方程的通解2. 求微分方程的通解3. 求微分方程组 在初始条件下的特解,并画出解函数的图形4. 分别用 ode23、ode45 求上述第 3 题中的微分方程初值问题的数值解(近似解),求解区间为利用画图来比较两种求解器之间的差异4 参考文献:1 Mastering Matlab 6, D. Hanselman, B. Littlefield, 清华大学出版,20022 赵静等编数学建模与数学实验(第3版)北京:高等教育出版社2008.3 姜启源编. 数学模型(第二版)北京:高等教育出版社.1993.4 石勇国. 蚂蚁追击问题与等角螺线. 宜宾学院学报. 2008,(6): 23-25.5 张伟年,杜正东,徐冰常微分方程北京:高等教育出版社20065 附录附录一:zhuji(x,y)的M文件function zhuji(x,y)clfv=1;dt=0.05;x(length(x)+1)=x(1);y(length(y)+1)=y(1);plot(x,y,'*') holdfor

温馨提示

  • 1. 本站所有资源如无特殊说明,都需要本地电脑安装OFFICE2007和PDF阅读器。图纸软件为CAD,CAXA,PROE,UG,SolidWorks等.压缩文件请下载最新的WinRAR软件解压。
  • 2. 本站的文档不包含任何第三方提供的附件图纸等,如果需要附件,请联系上传者。文件的所有权益归上传用户所有。
  • 3. 本站RAR压缩包中若带图纸,网页内容里面会有图纸预览,若没有图纸预览就没有图纸。
  • 4. 未经权益所有人同意不得将文件中的内容挪作商业或盈利用途。
  • 5. 人人文库网仅提供信息存储空间,仅对用户上传内容的表现方式做保护处理,对用户上传分享的文档内容本身不做任何修改或编辑,并不能对任何下载内容负责。
  • 6. 下载文件中如有侵权或不适当内容,请与我们联系,我们立即纠正。
  • 7. 本站不保证下载资源的准确性、安全性和完整性, 同时也不承担用户因使用这些下载资源对自己和他人造成任何形式的伤害或损失。

评论

0/150

提交评论