第五讲 Matlab求解微分方程

第五讲 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 是显式常微分方程：，在积分

⎪⎩

⎪

⎨

⎧

=

=

)

(

)

,(

y

t y

y

t

f

dt

dy

区间tspan=上，从到，用初始条件求解，

]

,

[

0f

t

t

t

f

t

y

要获得问题在其他指定时间点上的解，则令

,2

1

,

,t

t

t

tspan=（要求是单调的）．

]

,

,

,

[

,2

1

0f

t

t

t

t

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 方程；累计截断误差达 3)(x ∆大部分场合的首选算法

ode23 非刚性

单步算法；2、3阶Runge-Kutta

方程；累计截断误差达 3)(x ∆使用于精度较低的情形

ode113 非刚性 多步法；Adams 算法；高低精

度均可到

6

3

10~10

--计算时间比 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：求解微分方程

，并加以验证． 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)