卫星轨道常识.doc

5.14 常微分方程5.14.2 ODE解算指令的使用演示5.14.2.1 解算指令简洁格式使用示例【 * 例 5.14.2.1-1 】采用最简洁格式的 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 % 在主程序中定义全局变量传递参数 1G=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); %8X=YY(:,1); % 输出 Y 的第一列是位移数据 x(t)Y=YY(:,2); % 输出 Y 的第二列是位移数据 y(t)plot(X,Y,b,Linewidth,2); hold onaxis(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 卫星轨道 ? 【 * 例 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 flagcase % 按规定：这里必须是空串。在此为 真 时，完成以下导数计算。X=Y(1:2);V=Y(3:4);r=sqrt(sum(X.^2));Yd=[V;-G*ME*X/r^3];otherwiseerror([ Unknown flag flag . ]);end （2）按以下方法修改上例第（4）步中的程序，并运行之，可得相同结果。删去原主程序第 1 条指令；把原主程序的第 8 条指令改写为[t,YY]=ode45(DYDt2,tspan,Y0,[],G,ME); % 第 4 宗量取缺省设置% 第 5 、 6 宗量是被传递参数 5.14.2.2 解算指令较复杂格式的使用示例【 * 例 5.14.2.2-1 】带事件设置的 ODE 文件及主程序编写演示。本例将以较高精度计算卫星经过近地点和远地点的时间，并在图上标志。（1） ODE 文件的编写[DYDt3.m]function varargout=DYDt3(t,Y,flag,G,ME,tspan,Y0)% DYDt3.m 供主程序调用的 ODE 函数文件% 本文件自带三个子函数： f,fi,fev 。% t, Y 分别是自变量和一阶函数向量，是最基本的输入宗量。% flag 第三输入宗量，它专供解算指令（如 ode45 ）作调用通知。% 在运行中，解算指令会根据需要向 flag 发不同的字符串。% varargout 是 变 输出宗量。它由变维的元胞数组构成。每个元胞中可以存放指令% 所产生的任意形式的数据。switch flagcase