matlab常微分方程和常微分方程组的求解

下载的，感觉不错，共享一下

常微分方程和常微分方程组的求解

一、实验目的：

熟悉Matlab 软件中关于求解常微分方程和常微分方程组的各种命令，掌握利用Matlab 软件进行常微分方程和常微分方程组的求解。 二、相关知识

在MATLAB 中，由函数dsolve()解决常微分方程（组）的求解问题，其具体格式如下：

X=dsolve(‘eqn1’,’eqn2’,…)

函数dsolve 用来解符号常微分方程、方程组，如果没有初始条件，则求出通解，如果有初始条件，则求出特解。

例1：求解常微分方程1dy dx x y =

+的MATLAB 程序为：dsolve('Dy=1/(x+y)','x')，

注意，系统缺省的自变量为t ，因此这里要把自变量写明。 结果为：-lambertw(-C1*exp(-x-1))-x-1

其中：Y=lambertw(X)表示函数关系Y*exp(Y)=X 。

例2：求解常微分方程2

'''0yy y -=的MATLAB 程序为：

Y2=dsolve('y*D2y-Dy^2=0’,’x’) 结果为：

Y2 =[ exp((x+C2)/C1)]

[ C2]

我们看到有两个解，其中一个是常数。

例3：求常微分方程组253t

t dx x y e dt dy x y e dt ?++=???

?--=??通解的MATLAB 程序为： [X,Y]=dsolve('Dx+5*x+y=exp(t),Dy-x-3*y=exp(2*t)','t')

例4：求常微分方程组020

210cos ,224,0

t t t dx dy x t x dt dt dx dy y e y dt dt =-=?+-==???

?++==??通解的MATLAB 程序为：

[X,Y]=dsolve('Dx+2*x-Dy=10*cos(t),Dx+Dy+2*y=4*exp(-2*t)','x(0)=2','y(0)=0')

以上这些都是常微分方程的精确解法，也称为常微分方程的符号解。但是，我们知道，有大量的常微分方程虽然从理论上讲，其解是存在的，但我们却无法求出其解析解，此时，我们需要寻求方程的数值解，在求常微分方程数值解方面，MATLAB 具有丰富的函数，我们将其统称为solver ，其一般格式为：

[T,Y]=solver(odefun,tspan,y0)

该函数表示在区间tspan=[t0,tf]上，用初始条件y0求解显式常微分方程

'(,)y f t y =。

solver 为命令ode45，ode23，ode113，ode15s ，ode23s ，ode23t ，ode23tb 之一，这些命令各有特点。我们列表说明如下： 求解器 ODE 类型

特点

说明

ode45

非刚性

一步算法，4,5阶Runge-Kutta 方法累积截断误差3

()x ?

大部分场合的首选算

法

ode23 非刚性

一步算法，2,3阶Runge-Kutta 方法累积截断误差3()x ?

使用于精度较低的情

形

ode113

非刚性

多步法，Adams 算法，高低精

度均可达到3

6

10~10--

计算时间比ode45短 ode23t 适度刚性 采用梯形算法 适度刚性情形 ode15s 刚性 多步法，Gear’s 反向 数值积分，精度中等 若ode45失效时， 可尝试使用 ode23s

刚性

一步法，2阶Rosebrock 算法，

低精度。

当精度较低时， 计算时间比ode15s 短

odefun 为显式常微分方程'(,)y f t y =中的(,)f t y tspan 为求解区间，要获得问题在其他指定点012,,,

t t t 上的解，则令

012[,,,

,]f tspan t t t t =（要求i t 单调），

y0初始条件。

例5：求解常微分方程

2

'222y y x x =-++，00.5x ≤≤，(0)1y =的MATLAB 程序如下：fun=inline('-2*y+2*x*x+2*x')；[x,y]=ode23(fun,[0,0.5],1) 结果为：

x =

0,0.0400,0.0900,0.1400,0.1900,0.2400,0.2900,0.3400,0.3900,0.4400,0.4900,0.5000 y =

1.0000,0.9247,0.8434,0.7754,0.7199,0.6764,0.6440,0.6222,0.6105,0.6084,0.6154,0.6179

例6：求解常微分方程222

(1)0,(0)1,'(0)0d y dy y y y y dt dt μ--+===的解，并画

出解的图形。

分析：这是一个二阶非线性方程，用现成的方法均不能求解，但我们可以通过下面的变换，将二阶方程化为一阶方程组，即可求解。

令：1x y =，

2dy

x dt =

，7μ=，则得到：

接着，编写vdp.m 如下： function fy=vdp(t,x)

fy=[x(2);7*(1-x(1)^2)*x(2)-x(1)]; 再编写m 文件sy12_6.m 如下： y0=[1;0]

[t,x]=ode45(@vdp,[0,40],y0); y=x(:,1);dy=x(:,2); plot(t,y,t,dy) 三、实验内容

1．利用MATLAB 求常微分方程的初值问题38dy

y dx +=，02x y ==的解。

2．利用MATLAB 求常微分方程的初值问题2

(1)''2'x y xy +=，01x y ==，0'3

x y ==的解。

3．利用MATLAB 求常微分方程(4)

2'''''0y

y y -+=的解。

4．利用MATLAB 求常微分方程组003

24,2

30,0

t

t t dx dy x y e x dt dt dx x y y dt ==?++-==

???

?++==?

?的特解。

5．求解常微分方程2

''2(1)'0y y y y --+=，030x ≤≤，(0)1y =，'(0)0y =的特解，

并做出解函数的曲线图。