第五讲---Matlab求解微分方程教学提纲.doc

第五讲 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 是显式常微分方程：，在积分区间 tspan=上，从到，用初始条件求解，要获得问题在其他指定时间点上的解，则令 tspan=（要求是单调的）．
ezplot(x,y,[tmin,tmax])
符号函数的作图命令．x,y 为关于参数t 的符号函数，[tmin,tmax] 为 t 的取值范围．
inline()
建立一个内联函数．格式：inline('expr', 'var1', 'var2',…) ，注意括号里的表达式要加引号．
因为没有一种算法可以有效地解决所有的 ODE 问题，为此，Matlab 提供了多种求解器 Solver，对于不同的ODE 问题，采用不同的Solver．
求解器
Solver
ODE类型
特点
说明
ode45
非刚性
单步算法；4、5阶Runge-Kutta方程；累计截断误差达
大部分场合的首选算法
ode23
非刚性
单步算法；2、3阶Runge-Kutta方程；累计截断误差达
使用于精度较低的情形
ode113
非刚性
多步法；Adams算法；高低精度均可到
计算时间比 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：求解微分方程，并加以验证．
求解本问题的Matlab 程序为：
syms x y %line1
y=dsolve('Dy+2*x*y=x*exp(-x^2)','x') %line2
diff(y,x)+2 *x*y-x*exp(-x^2) %line3
simplify(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， 表明的确是微分方程的解．
例2：求微分方程在初始条件下的特解，并画出解函数的图形．
求解本问题的 Matlab 程序为：
syms x y
y=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.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.0000 0.0400 0.0900 0.1400 0.1900 0.2400
0.2900 0.3400 0.3900 0.4400 0.4900 0.5000
>> y'
ans =
1.0000 0.9247 0.8434 0.7754 0.7199 0.6764
0.6440 0.6222 0.6105 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,'*')

　　　图３导弹运行轨迹（解析法） 图４两种方法对比的导弹运行轨迹
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 on
x=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）建立平面直角坐标系．以正多边形的中心为原点, 设正多边形的一个顶点为起始点, 连接此点和原点作轴. 根据轴作出相应的轴; 选取足够小的进行采样．
（2）在每一时刻，计算每只蚂蚁在下一时刻时的坐标．不妨设甲追逐对象是乙，在时间时，甲的坐标为，乙的坐标为．甲在时在点(如图1), 其坐标为
，
其中. 同理，依次计算下一只蚂蚁在时的坐标．通过间隔进行采样，得到新一轮各个蚂蚁在一个新的正多边形位置坐标．
（4）重复2）步，直到充分小为止．
（5）连接每只蚂蚁在各时刻的位置，就得到所求的轨迹．

图1 在采样时间内，相连蚂蚁追击

用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. 求微分方程的通解．
2. 求微分方程的通解．
3. 求微分方程组

在初始条件下的特解，并画出解函数的图形．
4. 分别用 ode23、ode45 求上述第 3 题中的微分方程初值问题的数值解(近似解)，求解区间为．利用画图来比较两种求解器之间的差异．
4 参考文献:
[1] Mastering Matlab 6, D. Hanselman, B. Littlefield, 清华大学出版,2002
[2] 赵静等编．数学建模与数学实验（第3版）．北京：高等教育出版社．2008.
[3] 姜启源编. 数学模型（第二版）．北京：高等教育出版社.1993.
[4] 石勇国. 蚂蚁追击问题与等角螺线. 宜宾学院学报. 2008，(6): 23-25.
[5] 张伟年，杜正东，徐冰．常微分方程．北京：高等教育出版社．2006．
5 附录
附录一:zhuji(x,y)的M文件
function zhuji(x,y)
clf
v=1;
dt=0.05;
x(length(x)+1)=x(1);y(length(y)+1)=y(1);
plot(x,y,'*')
hold
for it=1:1000
for i=1:length(x)-1
d=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;
end
plot(x,y,'.')
hold on
x(length(x))=x(1); y(length(y))=y(1);
end