1、第五讲 Matlab 求解微分方程教学目的：学会用 MATLAB 求简单微分方程的解析解、数值解，加深对微分方程概念和应用的理解；针对一些具体的问题，如追击问题，掌握利用软件求解微分方程的过程；了解微分方程模型解决问题思维方法及技巧；体会微分方程建摸的艺术性教学重点：利用机理分析建模微分方程模型，掌握追击问题的建模方法，掌握利用 MATLAB 求解数值解教学难点：利用机理分析建模微分方程模型，通过举例，结合图形以及与恰当的假设突破教学难点1 微分方程相关函数（命令）及简介函数名 函数功能Dy 表示 y 关于自变量的一阶导数D2y 表示 y 关于自变量的二阶导数dsolve(equ1,equ2,
2、) 求微分方程的解析解，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 是显式常微分方程： ，0)(,ytfd在积分区间 tspan= 上，从 到 ，用初始条件,0f
3、t0tf求解，要获得问题在其他指定时间点 上的0y ,21t解，则令 tspan= （要求是单调的） ,210fttezplot(x,y,tmin,tmax) 符号函数的作图命令x,y 为关于参数 t 的符号函数，tmin,tmax 为 t 的取值范围inline() 建立一个内联函数格式：inline(expr, var1, var2,) ，注意括号里的表达式要加引号因为没有一种算法可以有效地解决所有的 ODE 问题，为此，Matlab 提供了多种求解器 Solver，对于不同的 ODE 问题，采用不同的 Solver求解器Solver ODE 类型 特点 说明ode45 非刚性单步算法；4
4、、5 阶 Runge-Kutta 方程；累计截断误差达3)(x大部分场合的首选算法ode23 非刚性单步算法；2、3 阶 Runge-Kutta 方程；累计截断误差达)(x使用于精度较低的情形ode113 非刚性多步法；Adams 算法；高低精度均可到 6310计算时间比 ode45 短ode23t 适度刚性 采用梯形算法 适度刚性情形ode15s 刚性 多步法；Gears 反向数值微分；精度中等 若 ode45 失效时，可尝试使用ode23s 刚性 单步法；2 阶 Rosebrock 算法；低精度 当精度较低时，计算时间比 ode15s 短ode23tb 刚性 梯形算法；低精度 当精度较低时
5、，计算时间比 ode15s 短要特别的是：ode23 、ode45 是极其常用的用来求解非刚性的标准形式的一阶常微分方程( 组) 的初值问题的解的 Matlab 的常用程序，其中：ode23 采用龙格-库塔 2 阶算法，用 3 阶公式作误差估计来调节步长，具有低等的精度ode45 则采用龙格-库塔 4 阶算法，用 5 阶公式作误差估计来调节步长，具有中等的精度2 求解微分方程的一些例子2.1 几个可以直接用 Matlab 求微分方程精确解的例子：例 1：求解微分方程 ，并加以验证2xeydx求解本问题的 Matlab 程序为：syms x y %line1y=dsolve(Dy+2*x*y=x
6、*exp(-x2),x) %line2diff(y,x)+2 *x*y-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(-x
7、2)*C1)(4) 行 line4 用 simplify() 函数对上式进行化简，结果为 0， 表明的确是微分方程的解)(xy例 2：求微分方程 在初始条件 下的特解，并画出解0xeyey2)1(函数的图形求解本问题的 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 格式)，即 ，xey此函数的图形如图 1：-6 -4 -2 0 2 4 6-30-20-1001020304050x1/x exp(x)+1/x exp(1
8、)图 1 y 关于 x 的函数图象2.2 用 ode23、ode45 等求解非刚性的标准形式的一阶常微分方程(组)的初值问题的数值解(近似解)例：求解微分方程初值问题 的数值解，求解范围为1)0(2yxdx区间0, 0.5fun=inline(-2*y+2*x2+2*x,x,y);x,y=ode23(fun,0,0.5,1); x; y;plot(x,y,o-) xans =0.0000 0.0400 0.0900 0.1400 0.1900 0.24000.2900 0.3400 0.3900 0.4400 0.4900 0.5000 yans =1.0000 0.9247 0.8434 0
9、.7754 0.7199 0.67640.6440 0.6222 0.6105 0.6084 0.6154 0.6179图形结果为图 20 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.50.60.650.70.750.80.850.90.951图 2 y 关于 x 的函数图像3 常微分在实际中的应用3.1 导弹追踪问题1、符号说明， ，乙舰的速率恒为 v0；设时刻 乙舰的坐标为 ，wt(),XtY导弹的坐标为 ；当零时刻，(),xty， 建立微分方程模型(0),0XY(,0)(,22 01(),10dykvdyxkxxw由微分方程模型解得 11()()
10、1,22()()kkxy k代入题设的数据 ，得到导弹的运行轨迹为/54655(1)()84yxx当 时 ，即当乙舰航行到点 处时被导弹击中. 被击中时间1x245y )245,1(为: . 若 v0=1, 则在 t=0.21 处被击中.0vt利用 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)()()(22yYxXdtyxX因乙舰以速度 沿直线 运动，设 ， ，因此导弹
11、0v101v5,1wt运动轨迹的参数方程为：0)(,)0( )(15)(22yxyttdt xxMATLAB 求解数值解程序如下，结果见图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)很明显数值计算的方法比较简单而且适用3.2 蚂蚁追击问题（1）建立平面直角坐标系以正多边形的中心为原点, 设正多边形的一个顶点为起始点, 连
12、接此点和原点作 轴. 根据 轴作出相应的 轴; 选取足够小xxy的 t进行采样（2）在每一时刻 ，计算每只蚂蚁在下一时刻 时的坐标不妨设甲追t t逐对象是乙，在时间 时，甲的坐标为 ，乙的坐标为 甲在A1(,)xyB2(,)xy时在 点(如图 1), 其坐标为tA，1(cossinxvtvt其中 . 同理，依次计算222211cos,sin,()()xydxyd下一只蚂蚁在 时的坐标通过间隔 t进行采样，得到新一轮各个蚂蚁在t一个新的正多边形位置坐标（4）重复 2）步，直到 充分小为止（5）连接每只蚂蚁在各时刻的位置，就得到所求的轨迹用 MALAB 求解并作图，函数 zhuJi(x,y)在附录
13、一定义， 如图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(2xyx2. 求微分方程 的通解ey53. 求微分方程组 A1,()xy2,()Bxyd图 1 在采样时间内，相连蚂蚁追击0yxdt在初始条件 下的特解，并画出解函数 的图形0|,1|0ttyx ()fx4. 分别用 ode23、ode45 求上述第 3 题中的微分方程初值问题的数值解(近似解 )，求解区间为 利用画图来比较两种求解器之间的差异,2t4 参考文献:1 Ma
14、stering 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 it=1:1000for i=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,.)hold onx(length(x)=x(1); y(length(y)=y(1);end
