第五讲 Matlab求解微分方程

第五讲 Matlab求解微分方程

教学目的：学会用MATLAB求简单微分方程的解析解、数值解，加深对微分方程概念和应用的理解；针对一些具体的问题，如追击问题，掌握利用软件求解微分方程的过程；了解微分方程模型解决问题思维方法及技巧；体会微分方程建摸的艺术性．

教学重点：利用机理分析建模微分方程模型，掌握追击问题的建模方法，掌握利用MATLAB求解数值解．

教学难点：利用机理分析建模微分方程模型，通过举例，结合图形以及与恰当的假设突破教学难点．

1微分方程相关函数（命令）及简介

因为没有一种算法可以有效地解决所有的ODE 问题，为此，Matlab 提供了多种求解器Solver，对于不同的ODE 问题，采用不同的Solver．

阶常微分方程(组)的初值问题的解的 Matlab 的常用程序，其中：

ode23 采用龙格-库塔2 阶算法，用3 阶公式作误差估计来调节步长，具有低等的精度．

ode45 则采用龙格-库塔4 阶算法，用5 阶公式作误差估计来调节步长，具有中等的精度．

2 求解微分方程的一些例子

几个可以直接用 Matlab 求微分方程精确解的例子：

例1：求解微分方程

22x xe xy dx

dy

-=+，并加以验证． 求解本问题的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， 表明)(x y y =的确是微分方程的解．

例2：求微分方程0'=-+x e y xy 在初始条件e y 2)1(=下的特解，并画出解函数的图形．

求解本问题的 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 格式)，即x

e e y x

+=，

此函数的图形如图 1：

图1 y 关于x 的函数图象

用ode23、ode45等求解非刚性的标准形式的一阶常微分方程(组)的初值问题的数值解(近似解)．

例３：求解微分方程初值问题⎪⎩⎪⎨⎧=++-=1)0(2222

y x

x y dx dy 的数值解，求解范围为

区间[0, ]．

fun=inline('-2*y+2*x^2+2*x','x','y'); [x,y]=ode23(fun,[0,],1); x'; y'; plot(x,y,'o-') >> x' ans =

>> y' ans =

图形结果为图 2．

图2 y 关于x 的函数图像

3 常微分在实际中的应用 导弹追踪问题

1、符号说明，w ，乙舰的速率恒为v 0；设时刻t 乙舰的坐标为((),())X t Y t ，导弹的坐标为((),())x t y t ；当零时刻，((0),(0))(1,0)X Y =，((0),(0))(0,0)x y =．建立微分方程模型．

20

2

,01(0)0,'(0)0v d y

k x dx w y y ⎧⎪⎪= =<<⎨⎪⎪==⎩

由微分方程模型解得

11(1)(1)11

(),12(1)2(1)2(1)2(1)

k k x x y x k k k k k +-+--=-+-≠+-+-++

代入题设的数据1/5k =，得到导弹的运行轨迹为

46

55555(1)(1)81224

y x x =--+-+

当1=x 时245=

y ，即当乙舰航行到点)24

5

,1(处时被导弹击中. 被击中时间

为:0

245

v v y t ==

. 若v 0=1, 则在t =处被击中. 利用MALAB 作图如图３．

clear, x=0::1; y=-5*(1-x).^(4/5)/8+5*(1-x).^(6/5)/12+5/24; plot(x,y,'*')

图３导弹运行轨迹（解析法） 图４两种方法对比的导弹运行轨迹

2、数值方法求解．设导弹速率恒为w ，则得到参数方程为

⎪⎪⎩

⎪

⎪⎨⎧--+-=--+-=)

()()()()()(222

2y Y y Y x X w

dt dy x X y Y x X w dt dx

因乙舰以速度0v 沿直线1x =运动，设01v =，5,1,w X Y t ===，因此导弹运动轨迹的参数方程为：

⎪⎪⎪⎩

⎪⎪

⎪⎨

⎧==--+-=--+-=0)0(,0)0()()()1(5)1()()1(52

22

2y x y t y t x dt

dy

x y t x dt dx

MATLAB 求解数值解程序如下，结果见图４ t0=0,tf=;

[t,y]=ode45('eq2',[t0 tf],[0 0]); X=1;Y=0::;plot(X,Y,'-')

plot(y(:,1),y(:,2),'*'),hold on

x=0::1; y=-5*(1-x).^(4/5)/8+5*(1-x).^(6/5)/12+5/24; plot(x,y,'r')

很明显数值计算的方法比较简单而且适用． 蚂蚁追击问题

（1）建立平面直角坐标系．以正多边形的中心为原点, 设正多边形的一个顶点为起始点, 连接此点和原点作x 轴. 根据x 轴作出相应的y 轴; 选取足够小的进行采样．

（2）在每一时刻t ，计算每只蚂蚁在下一时刻t t +∆时的坐标．不妨设甲追