matlab用dsolv求解常微分方程

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

r = dsolve('eq1,eq2,...', 'cond1,cond2,...', 'v')

'eq1,eq2,...'为微分方程或微分方程组，'cond1,cond2,...',是初始条件或边界条件，'v'是变量，默认的变量是't'。

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

例1：求解常微分方程的MATLAB程序为：dsolve('Dy=1/(x+y)','x')  ，

注意，系统缺省的自变量为t，因此这里要把自变量写明。

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

例2：求解常微分方程的MATLAB程序为：

Y2=dsolve('y*D2y-Dy^2=0','x')  

Y2=dsolve('D2y*y-Dy^2=0','x')  

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

例3：求常微分方程组通解的MATLAB程序为：

[X,Y]=dsolve('Dx+5*x+y=exp(t),Dy-x-3*y=exp(2*t)','t')  

例4：求常微分方程组通解的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求解显式常微分方程。

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

tspan为求解区间，要获得问题在其他指定点上的解，则令（要求单调递增或递减），y0初始条件。

例5：求解常微分方程，，的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：求解常微分方程的解，并画出解的图形。

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

令：，，，则得到：

解：

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）问题的形成

轨道上的卫星，在牛顿第二定律，和万有引力定律作用下有

，引力常数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

（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.e7;                        %

XE=RE*XE;YE=RE*YE;ZE=0*ZE;        %

mesh(XE,YE,ZE),hold off            %  

练习：

1．利用MATLAB求常微分方程的初值问题，的解。

r=dsolve('Dy+3*y=0','y(0)=2','x')   

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

r=dsolve('D2y*(1+x^2)-2*x*Dy=0','y(0)=1,Dy(0)=3','x')  

3．利用MATLAB求常微分方程的解。

解：y=dsolve('D4y-2*D3y+D2y','x')  

4．利用MATLAB求常微分方程组的特解。

[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．求解常微分方程，，，的特解，并作出解函数的曲线图。

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)

legend('y','dy')  

求解特殊函数方程

勒让德方程的求解

r=dsolve('(1-x^2)*D2y-2*x*Dy+y*l*(l+1)=0','x')  

连带勒让德方程的求解：

r=dsolve('(1-x^2)*D2y-2*x*Dy+y*(l*(l+1)-m^2/(1-x^2))=0','x')  

贝塞尔方程

r=dsolve('x^2*D2y+x*Dy+(x^2-v^2)*y=0','x')  

利用maple

maple dsolve((1-x^2)*diff(y(x),x$2)-2*x*diff(y(x),x)+n*(n+1)*y(x) = 0, y(x))  

利用MAPLE的深层符号计算资源

经典特殊函数的调用

MAPLE库函数在线帮助的检索树

发挥MAPLE的计算潜力

调用MAPLE函数

【例5.7.3.1-1】求递推方程的通解。

（1）

gs1=maple('rsolve(f(n)=-3*f(n-1)-2*f(n-2),f(k));')  

（2）调用格式二

gs2=maple('rsolve','f(n)=-3*f(n-1)-2*f(n-2)','f(k)')  

【例5.7.3.1-2】求的Hessian矩阵。

（1）

FH1=maple('hessian(x*y*z,[x,y,z]);')  

（2）

FH2=maple('hessian','x*y*z','[x,y,z]')  

（3）

FH=sym(FH2)  

【例5.7.3.1-3】求在处展开的截断8阶小量的台劳近似式。

（1）

TL1=maple('mtaylor(sin(x^2+y^2),[x=0,y=0],8)')  

（2）

maple('readlib(mtaylor);');

TL2=maple('mtaylor(sin(x^2+y^2),[x=0,y=0],8)');

pretty(sym(TL2))  

运行MAPLE程序

【例5.7.3.2-1】目标：设计求取一般隐函数的导数解析解的程序，并要求：该程序能象MAPLE原有函数一样可被永久调用。

（1）

[DYDZZY.src]

DYDZZY:=proc(f)

# DYDZZY(f) is used to get the derivate of

#         an implicit function

local Eq,deq,imderiv;

Eq:='Eq';

Eq:=f;

deq:=diff(Eq,x);

readlib(isolate);

imderiv:=isolate(deq,diff(y(x),x));

end;

（2）

procread('DYDZZY.src')  

（3）

s1=maple('DYDZZY(x=log(x+y(x)));')

s2=maple('DYDZZY(x^2*y(x)-exp(2*x)=sin(y(x)))')

s3=maple('DYDZZY','cos(x+sin(y(x)))=sin(y(x))')  

（4）

clear maplemex

procread('DYDZZY.src');    

maple('save(`DYDZZY.m`)');  

（5）

maple('read','`DYDZZY.m`');

ss2=maple('DYDZZY(x^2*y(x)-exp(2*x)=sin(y(x)))')