matlab求解常微分方程.docx

用matlab求解常微分方程

在MATLAB中，由函数dsolve()解决常微分方程(组)的求解问题，其具体格式如下：r 二dsolve('eql,eq2,・••字condl,cond2,・.・；V)

匕ql,eq2,・・・*为微分方程或微分方程组，，condl,cond2,.・・；是初始条件或边界条件，P是独立变量，默认的独立变量是讥

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

dy _1

例1：求解常微分方程莎一石的MATLAB程序为：

dsolve(* Dy=l/(x+y) 1r!x1),

注意，系统缺省的自变量为t,因此这里要把自变量写明。其中：Y=lambertw(X)表示函数关系Y*exp(Y)二X。

例2：求解常微分方程E'-y— 0的MATLAB程序为：

Y2=dsolve(1y*D2y-Dy A2=01, 1x f)

Y2=dsolve(!D2y*y-Dy A2=0 J )

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

dx 心 ？ —+ 5x + y = e dt

空_兀_3『= g2f 例3：求常微分方程组〔力 ' 通解的MATLAB 程序为： [X,Y]=dsolve(f Dx+5*x+y=exp(t),Dy-x-3*y=exp(2*t) 1, 111) [X,Y]=dsolve(1Dx+2 *x-Dy=l0 * cos(t),Dx+Dy+2 *y=4 *exp(-

2*t) T ,f x(0)=2,y(0)=0f ,f t T

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

i°cosr, 7

=2 例4：求常微分方程组 y = 0 z 通解的MATLAB 程序为:

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

该函数表示在区间tspan=[tO,tf]±,用初始条件yO求解显式常微分方程卩="，刃。

solver 为命令odc45, odc23, ode 113, ode 15s, ode23s, ode23t, ode23tb 之一,这些命令各有特点。我们列表说明如下:

odefun为显式常微分方程）,”，刃中的“丿）

tspan为求解区间，要获得问题在其他指定点心，也，…上的解，则令切期二［心,/“，…心］（要求-单调递增或递减），yO初始条件。

例5：求解常微分方程y' = -2y + 2x2+2x, 0

y=dsolve（1Dy=-2*y+2*x A2+2*x1, 1y （0）=11, 1x!） x=0:0.01:0.5;

yy=subs(y,x);

fun=inline (*-

2*y+2*x*x+2*x1); [x,y]=odel5s(fun, [0:0.01:0.5],1);ys=x ・ *x+exp(一

2*x);

Plot(x,y, f r f,x,ys, f b f )

例6：求解常微分方程分冷+尸°"3)"'八°2()的解，并画出解的图形。

分析：这是一个二阶非线性方程(函数以及所有偏导数军委一次幕的是现性方程，高于一次的为非线性方程)，用现成的方法均不能求解，但我们可以通过下面的变换，将二阶方程化为一阶方程组，即可求解。

= dy_

令：兀产儿氐-亦，“ =7,则得到：

~~ =兀2 ，兀| (0) = 1

at

<

= 7(1-%,2 )x2 - Xj, x2 (0) = 0

dt

解:

function [dfy]=mytt(t,fy)

%fl=y;f2=dy/dt

%求二阶非线性微分方程时，把一阶、二阶直到(ml)阶导数用另外一个函数代替

%用ode45命令时，必须表示成Y』f(t,Y)的形式

%Y=[yl;y2;y3],Y,=[yl\y2\y3>[y2;y3;f(yl,y2,y3)],

%其中yl=y,y2=y',y3=y"

%更高阶时类似

dfy=[fy(2);7*(l ・fy ⑴八2)*fy(2)・fy ⑴]；

clear;clc

[t,yy]=ode45Cmytt\[0 40],[1;0]);

plot(t,yy)

legend('yTdy')

【例4.14.2.1-1］采用ODE解算扌旨令研究围绕地球旋转的卫星轨道。

(1)问题的形成

轨道上的卫星，在牛顿第二定律F=ma = m^,和万有引力定律F T叫；作用下有dt r

2

a^L = -G^r ,引力常数G=6.672*10 H(N.m2/kg2) ,M

E

=5.97*1024(kg)是地球的质量。假dr r

定卫星以初速度v y(0)=4000m/s在x(0)=4.2*l()7(m)处进入轨道。

(2) 构成一阶微分方程组

令 Y 二[刃力力为]「二[尤 V V x v y ]T =[x y * y]T

儿 (…严 儿 (宀)计

(3) 根据上式

[dYdt.m]

function Ycl=DYdt (t, Y)

% t

% Y

global G ME

% xy 二Y(l:2) ;Vxy=Y(3:4) ; %

r=sqrt(sum(xy.八

2)); Yd 二[Vxy ;-G*ME*xy//3] ; %

(4)

global G ME % <1>

G=6 ・ 672e-ll;ME=5 ・ 97e24;vy0=4000;x0=-4 ・ 2e7;t0=0;tf=60*60*24*9; tspan=[tO,tf];

%

-GM E —G

M