matlab解无解析解微分方程组,[转载]Matlab解常微分方程深入解析

采用最简洁格式的ODE文件和解算指令，研究围绕地球旋转的卫星轨道。

图 5.14.2.1-1 地球轨道卫星运动加速度

(1)问题的形成

轨道上运动的卫星，在Newton第二定律 ，和万有引力定律 作用下，有 。即 ， ，而 。这里 是引力常数，

是地球的质量。又假定卫星以初速度 在 处进入轨道。

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

令 ，则

(5.14.2.1-5)

初始条件为

(5.14.2.1-6)

(3)根据式(5.14.2.1-5)编写最简洁格式的ODE文件

[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]; %Yd必须按式(5.14.2.1-5)编写，是与Y同维的列向量。

(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]; %按式(5.14.2.1-6)给定初值向量

[t,YY]=ode45('DYDt',tspan,Y0); %<8>

X=YY(:,1); %输出Y的第一列是位移数据x(t)

Y=YY(:,2); %输出Y的第二列是位移数据y(t)

plot(X,Y,'b','Linewidth',2); hold on

axis('image') %保证x、y轴等长刻度，且坐标框恰包容图形

[XE,YE,ZE] =

sphere(10); %产生单位球面数据(三维坐标)

RE=0.64e7; %地球半径

XE=RE*XE;YE=RE*YE;ZE=0*ZE; %坐标纸上的地球平面数据

mesh(XE,YE,ZE),hold

off %绘地球示意图

图 5.14.2.1-2

卫星轨道

上例中，程序间的参数(如G和ME)传送，是依靠全局变量形式实现的。一般说来，编写程序时，应尽量少用全局变量，以免引起混乱。本例演示参数如何在指令间直接传送。

要实现参数直接传送，必须对上例中的程序进行修改，具体如下：

(1)重写ODE文件

[DYDt2.m]

function Yd=DYDt2(t,Y,flag,G,ME)

% flag 按ODE文件格式规定，必须是第三输入宗量。对它的赋值由ode45指令自动产生。

% 第4、5宗量是被传递的参数

switch flag

case

'' %按规定：这里必须是空串。在此为"真"时，完成以下导数计算。

X=Y(1:2);V=Y(3:4);r=sqrt(sum(X.^2));Yd=[V;-G*ME*X/r^3];

otherwise

end

(2)按以下方法修改上例第(4)步中的程序，并运行之，可得相同结果。

删去原主程序第<1>条指令；把原主程序的第<8>条指令改写为

[t,YY]=ode45('DYDt2',tspan,Y0,[],G,ME); %第4宗量取缺省设置

%第5、6宗量是被传递参数

解算指令较复杂格式的使用示例

带事件设置的ODE文件及主程序编写演示。本例将以较高精度计算卫星经过近地点和远地点的时间，并在图上标志。

(1)ODE文件的编写

[DYDt3.m]

function varargout=DYDt3(t,