用MATLAB解常微分方程

实验四 求微分方程的解

一、问题背景与实验目的

实际应用问题通过数学建模所归纳而得到的方程,绝大多数都是微分方

程,真正能得到代数方程的机会很少.另一方面,能够求解的微分方程也是十分有限的,特别是高阶方程和偏微分方程(组).这就要求我们必须研究微分方程(组)精品文档,你值得期待

的解法,既要研究微分方程(组)的解析解法(精确解),更要研究微分方程(组)的数值解法(近似解).

对微分方程(组)的解析解法(精确解),Matlab 有专门的函数可以用,本实验将作一定的介绍.

本实验将主要研究微分方程(组)的数值解法(近似解),重点介绍 Euler 折线法.

二、相关函数(命令)及简介

1.dsolve('equ1','equ2',…):Matlab 求微分方程的解析解.equ1、equ2、…为方程(或条件).写方程(或条件)时用 Dy 表示y 关于自变量的一阶导数,用用 D2y 表示 y 关于自变量的二阶导数,依此类推.

2.simplify(s):对表达式 s 使用 maple 的化简规则进行化简. 例如: syms x

simplify(sin(x)^2 + cos(x)^2) ans=1

3.[r,how]=simple(s):由于 Matlab 提供了多种化简规则,simple 命令就是对表达式 s 用各种规则进行化简,然后用 r 返回最简形式,how 返回形成这种形式所用的规则.

例如: syms x

[r,how]=simple(cos(x)^2-sin(x)^2) r = cos(2*x) how = combine

4.[T,Y] = solver(odefun,tspan,y0) 求微分方程的数值解. 说明:

(1) 其中的 solver为命令 ode45、ode23、ode113、ode15s、ode23s、ode23t、ode23tb 之一.

?dy??f(t,y)(2) odefun 是显式常微分方程:?dt

??y(t0)?y0(3) 在积分区间 tspan=[t0,tf]上,从t0到tf,用初始条件y0求解.

(4) 要获得问题在其他指定时间点t0,t1,t2,?上的解,则令 tspan=

[t0,t1,t2,?,tf](要求是单调的).

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

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

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

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

5.ezplot(x,y,[tmin,tmax]):符号函数的作图命令.x,y 为关于参数t 的符号函数,[tmin,tmax] 为 t 的取值范围.

6.inline():建立一个内联函数.格式:inline('expr', 'var1', 'var2',…) ,注意括号里的表达式要加引号.

例:Q = dblquad(inline('y*sin(x)'), pi, 2*pi, 0, pi)

三、实验内容

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

2dy例1:求解微分方程?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)

(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:

联系客服:779662525#qq.com(#替换为@) 苏ICP备20003344号-4