第五讲 Matlab求解微分方程

第五讲 Matlab求解微分方程

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

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

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

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

函数功能 Dy 表示y 关于自变量的一阶导数 D2y 表示 y 关于自变量的二阶导数 dsolve('equ1','equ2',…) 求微分方程的解析解，equ1、equ2、…为方程（或条件） simplify(s) 对表达式 s 使用 maple 的化简规则进行化简 simple 命令就是对表达式 s 用各种规则进行化简，然后[r,how]=simple(s) 用 r 返回最简形式，how 返回形成这种形式所用的规则． 求微分方程的数值解，其中的 solver为命令 ode45、 ode23、ode113、ode15s、ode23s、ode23t、ode23tb 之一， ?dy??f(t,y)[T,Y] = odefun 是显式常微分方程：?dt，在积分区间 ?solver(odefun,tspan,y0) ?y(t0)?y0tspan=[t0,tf]上，从t0到tf，用初始条件y0求解，要获得问题在其他指定时间点t0,t1,t2,?上的解，则令 tspan=[t0,t1,t2,?,tf]（要求是单调的）． ezplot(x,y,[tmin,tmax]) 符号函数的作图命令．x,y 为关于参数t 的符号函数，[tmin,tmax] 为 t 的取值范围． 建立一个内联函数．格式：inline('expr', 'var1', 'var2',…) ，inline() 注意括号里的表达式要加引号． 函数名 因为没有一种算法可以有效地解决所有的 ODE 问题，为此，Matlab 提供了多种求解器 Solver，对于不同的ODE 问题，采用不同的Solver．

求解器 Solver ode45 ODE类型 非刚性 特点 单步算法；4、5阶Runge-Kutta方程；累计截断误差达(?x)3 单步算法；2、3阶Runge-Kutta说明 大部分场合的首选算法 ode23 非刚性 方程；累计截断误差达(?x)3 多步法；Adams算法；高低精使用于精度较低的情形 ode113 ode23t 非刚性 适度刚性 度均可到10?3~10 ?6计算时间比 ode45 短 适度刚性情形 若 ode45 失效时，可尝ode15s 刚性 试使用 当精度较低时，计算时ode23s 刚性 间比 ode15s 短 当精度较低时，计算时ode23tb 刚性 梯形算法；低精度 间比 ode15s 短 要特别的是：ode23、ode45 是极其常用的用来求解非刚性的标准形式的一阶常微分方程(组)的初值问题的解的 Matlab 的常用程序，其中：

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

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

采用梯形算法 多步法；Gear's反向数值微分；精度中等 单步法；2阶 Rosebrock 算法；低精度 2 求解微分方程的一些例子

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

2dy?2xy?xe?x，并加以验证． dx求解本问题的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)

例1：求解微分方程

(4) 行line4 用 simplify() 函数对上式进行化简，结果为 0， 表明y?y(x)的确是微分方程的解．

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

求解本问题的 Matlab 程序为： syms x y

y=dsolve('x*Dy+y-exp(x)=0','y(1)=2*exp(1)','x') ezplot(y)

e?ex微分方程的特解为：y=1/x*exp(x)+1/x* exp (1) (Matlab格式)，即y?，

x此函数的图形如图 1：

1/x exp(x)+1/x exp(1)50403020100-10-20-30-6-4-20x246 图1 y关于x的函数图象

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

?dy2???2y?2x?2x例３：求解微分方程初值问题?dx的数值解，求解范围为

??y(0)?1区间[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．

10.950.90.850.80.750.70.650.600.050.10.150.20.250.30.350.40.450.5 图2 y关于x的函数图像

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

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

?dy2k1?()?d2y?dx,????k?v0?,0?x?1 ?2dx(1?x)w??y(0)?0,y'(0)?0?由微分方程模型解得

(1?x)k?1(1?x)?k?111y(x)????,k?1

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

46555y??(1?x)5?(1?x)5?

81224当x?1时y?55，即当乙舰航行到点(1, )处时被导弹击中. 被击中时间2424