matlab求解常微分方程
- 格式：pdf
- 大小：226.07 KB
- 文档页数：15
用matlab求解常微分方程
在MATLAB中，由函数dsolve()解决常微分方程（组）的求解问题，其具体格式如
下：
r = dsolve('eq1,eq2,...', 'cond1,cond2,...', 'v')
'eq1,eq2,...'为微分方程或微分方程组，'cond1,cond2,...',是初始条件或边界条件，'v'是
独立变量，默认的独立变量是't'。
函数dsolve用来解符号常微分方程、方程组，如果没有初始条件，则求出通解，如
果有初始条件，则求出特解。
例1：求解常微分方程1dydxxy=+的MATLAB程序为：dsolve('Dy=1/(x+y)','x') ，
注意，系统缺省的自变量为t，因此这里要把自变量写明。
其中：Y=lambertw(X)表示函数关系Y*exp(Y)=X。
例2：求解常微分方程的MATLAB程序为： 2'''0yyy−=Y2=dsolve('y*D2y-Dy^2=0','x')
Y2=dsolve('D2y*y-Dy^2=0','x')
我们看到有两个解，其中一个是常数0。
例3
：求常微分方程组25
3t
tdxxyedtdyxyedt⎧++=⎪⎪⎨⎪−−=⎪⎩通解的MATLAB程序为： [X,Y]=dsolve('Dx+5*x+y=exp(t),Dy-x-3*y=exp(2*t)','t')
例4
：求常微分方程组020210cos,
24,tttdxdyxtxdtdtdxdyyeydtdt=−=⎧+−==⎪⎪⎨⎪++==⎪⎩2
0通解的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','t')
以上这些都是常微分方程的精确解法，也称为常微分方程的符号解。但是，我们知
道，有大量的常微分方程虽然从理论上讲，其解是存在的，但我们却无法求出其解析
解，此时，我们需要寻求方程的数值解，在求常微分方程数值解方面，MATLAB具有丰
富的函数，我们将其统称为solver，其一般格式为： [T,Y]=solver(odefun,tspan,y0)
该函数表示在区间tspan=[t0,tf]上，用初始条件y0求解显式常微分方程。 '(,yfty=)
solver为命令ode45，ode23，ode113，ode15s，ode23s，ode23t，ode23tb之一，这些
命令各有特点。我们列表说明如下： 求解器 特点 说明
ode45 一步算法，4,5阶Runge-Kutta 方法累积截断误差3()xΔ大部分场合的首选算法
ode23 一步算法，2,3阶Runge-Kutta 方法累积截断误差3()xΔ使用于精度较低的情形
ode113多步法，Adams算法，高低精度均可达到3610~10−−计算时间比ode45短
ode23t 采用梯形算法 适度刚性情形
ode15s 多步法，Gear’s反向 数值积分，精度中等 若ode45失效时， 可尝试使用
ode23s 一步法，2阶Rosebrock算法， 低精度。 当精度较低时， 计算时间比ode15s短 odefun为显式常微分方程中的'(,yfty=)(,)fty
tspan为求解区间，要获得问题在其他指定点上的解，则令012,,,ttt012[,,,,]ftspantttt=
（要求it单调递增或递减），y0初始条件。
例5：求解常微分方程2'222yyx=−++x.5，00x≤≤，(0)1y=的MATLAB程序如下： y=dsolve('Dy=-2*y+2*x^2+2*x','y(0)=1','x')
x=0:0.01:0.5; yy=subs(y,x);
fun=inline('-
2*y+2*x*x+2*x');[x,y]=ode15s(fun,[0:0.01:0.5],1);ys=x.*x+exp(-2*x);
plot(x,y,'r',x,ys,'b')
例6：求解常微分方程222(1)0,(0)1,'(0)0dydyyyyydtdtμ−−+===的解，并画出解的图形。 分析：这是一个二阶非线性方程(函数以及所有偏导数军委一次幂的是现性方程，高
于一次的为非线性方程)，用现成的方法均不能求解，但我们可以通过下面的变换，将二
阶方程化为一阶方程组，即可求解。
令：1xy=，2dyxdt=，7μ=，则得到：
121221212,(0)1
7(1),(0)0dxxxdtdxxxxxdt⎧==⎪⎪⎨⎪=−−=⎪⎩ 解： function [dfy]=mytt(t,fy)
%f1=y;f2=dy/dt
%求二阶非线性微分方程时，把一阶、二阶直到(n-1)阶导数用另外一个函数代替
%用ode45命令时，必须表示成Y'=f(t,Y)的形式
%Y=[y1;y2;y3],Y'=[y1';y2';y3']=[y2;y3;f(y1,y2,y3)],
%其中y1=y,y2=y',y3=y''
%更高阶时类似
dfy=[fy(2);7*(1-fy(1)^2)*fy(2)-fy(1)];
clear;clc
[t,yy]=ode45('mytt',[0 40],[1;0]);
plot(t,yy)
legend('y','dy')
【例4.14.2.1-1】采用ODE解算指令研究围绕地球旋转的卫星轨道。
（1）问题的形成 轨道上的卫星，在牛顿第二定律22drFmamdt==，和万有引力定律3EmMFGr−=−r作用下有
22EMdraGdtr−=−3r ，引力常数G=6.672*10-11(N.m2/kg2) ,ME=5.97*1024(kg)是地球的质量。假
定卫星以初速度vy(0)=4000m/s在x(0)=-4.2*107(m)处进入轨道。 （2）构成一阶微分方程组
令Y=[y1 y2 y3 y4]T=[x y vx vy]T=[x y x' y']T
31421223/2342223/2'''()'()'()xyExyEyvyyxyyYtGMayxyayyGMxy⎡⎤⎢⎥⎡⎤⎡⎤⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥===−⎢⎥⎢⎥⎢⎥+⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦⎣⎦−⎢⎥+⎣⎦i
i
（3）根据上式
[dYdt.m]
function Yd=DYdt(t,Y)
% t
% Y
global G ME %
xy=Y(1:2);Vxy=Y(3:4); %
r=sqrt(sum(xy.^2));
Yd=[Vxy;-G*ME*xy/r^3]; %
（4） global G ME % <1>
G=6.672e-11;ME=5.97e24;vy0=4000;x0=-4.2e7;t0=0;tf=60*60*24*9;
tspan=[t0,tf]; %
Y0=[x0;0;0;vy0]; %
[t,YY]=ode45('DYdt',tspan,Y0);% <8> X=YY(:,1); %
Y=YY(:,2); %
plot(X,Y,'b','Linewidth',2); hold on
%axis('image') %
[XE,YE,ZE] = sphere(10); %
RE=0.64e7; %
XE=RE*XE;YE=RE*YE;ZE=0*ZE; %
mesh(XE,YE,ZE),hold off %
练习：
1．利用MATLAB求常微分方程的初值问题38dyydx+=，02xy==的解。 r=dsolve('Dy+3*y=8','y(0)=2','x')
function dy=myddy(t,Y)
%
dy=[-3*Y(1)+8];
[t,yy]=ode45('myddy',[0:0.01:10],[2]); yys=subs(r,t);
plot(t,yy,t,yys);
legend('y','yys')
2．利用MATLAB求常微分方程的初值问题2(1)''2'xyx+=y，01xy==，0'xy==3的解。 r=dsolve('D2y*(1+x^2)-2*x*Dy=0','y(0)=1,Dy(0)=3','x')
%%% function yy=mydy2(x,Y) % yy=[Y(2);Y(2)*2*x/(1+x^2)];
clear;clc [t,YY]=ode45('mydy2',[0 30],[1;3]); ys=1+t.*t.*t+3*t; plot(t,YY,t,ys) legend('y','dy','ys') 3．利用MATLAB求常微分方程的解。 (4)2'''''0yyy−+=解：y=dsolve('D4y-2*D3y+D2y','x')
4．利用MATLAB
求常微分方程组0
0324,230,tt
tdxdyxyexdtdtdxxyydt=
=⎧
0++−==⎪⎪⎨⎪++==⎪⎩的特解。 [X,Y]=dsolve('Dx*2+4*x+Dy-
y=exp(t),Dx+3*x+y=0','x(0)=1.5,y(0)=0','t')
5．求解常微分方程，2''2(1)'0yyyy−−+=030x≤≤，(0)1y=，'(0)0y=的特解，并作出解函数的曲线图。 r=dsolve('D2y-2*(1-y^2)*Dy+y=0','y(0)=0,Dy(0)=0','x')
function DY=mytt2(t,Y)
DY=[Y(2);2*(1-Y(1)^2)*Y(2)+Y(1)];
clear;clc
[t,YY]=ode45('mytt2',[0 30],[1;0]);
plot(t,YY)
