学号：姓名：
《Matlab/Simulink在数学计算与仿真中的应用》大作业1.假设地球和火星绕太阳运转的半径分别为r和2r，利用comet指令动画显示从地球
到火星的转移轨迹（r可以任意取值，要求实时显示探测器、太阳、地球和火星的位置）。
解函数function comet(varargin)
[ax,args,nargs] = axescheck(varargin{:});
error(nargchk(1,3,nargs,'struct'));
% Parse the rest of the inputs
if nargs < 2, x = args{1}; y = x; x = 1:length(y); end
if nargs == 2, [x,y] = deal(args{:}); end
if nargs < 3, p = 0.10; end
if nargs == 3, [x,y,p] = deal(args{:}); end
if ~isscalar(p) || ~isreal(p) || p < 0 || p >= 1
error('MATLAB:comet:InvalidP', ...
'The input ''p'' must be a real scalar between 0 and 1.'); End
指令 %particle_motion
t = 0:5:16013;
r1=6.7e6;%随便给定参数
%---------------------------
r2=2*r1;
g=9.8;
R=6.378e6;
m=g*R^2;
%内轨道
v_inner=sqrt(m/r1);
w_inner=v_inner/r1;
x_inter=r1*cos(w_inner*t);
y_inter=r1*sin(w_inner*t);
%外轨道
v_outer=sqrt(m/r2);
w_outer=v_outer/r2;
x_outer=r2*cos(w_outer*t);
y_outer=r2*sin(w_outer*t);
%控制器转移轨道
a=(r1+r2)/2;
E=-m/(2*a);
V_near=sqrt(m*(2/r1-2/(r1+r2)));%转移轨道在近地点的速度
V_far=sqrt(m*(2/r2-2/(r1+r2)));%转移轨道在远地点的速度
h=r1*V_near;%由于在近地点的速度垂直于位置失量, h是转移轨道的比动量矩
e=sqrt(1+2*E*h^2/m^2);%e为椭圆轨迹的偏心率
TOF=pi*sqrt(a^3/m);%转移轨道是椭圆的一半及飞行时间是周期的一半(开普勒第三定律)
w=pi/TOF;%椭圆轨迹的角速度
c=a*e;
b=sqrt(a^2-c^2);
x_ellipse=a*cos(w*t)-0.5*r1;
y_ellipse=b*sin(w*t);
%动画显示运动轨迹
x=[x_inter;x_outer;x_ellipse]';
y=[y_inter;y_outer;y_ellipse]';
comet(x,y)
%---------------------------
动态图像如下：
2．利用两种不同途径求解边值问题df
dx
f g
dg
dx
f g f g
=+=-+==
34430001
,,(),()的
解.途径一：指令syms f g
[f,g]=dsolve('Df=3*f+4*g,Dg=-4*f+3*g','f(0)=0,g(0)=1');
disp('f=');disp(f)
disp('g=');disp(g)
结果（Matlab 7.8版本）f=
i/(2*exp(t*(4*i - 3))) - (i*exp(t*(4*i + 3)))/2
g=exp(t*(4*i + 3))/2 + 1/(2*exp(t*(4*i - 3)))
（Matlab 6.5版本）
f=exp(3*t)*sin(4*t)
g=exp(3*t)*cos(4*t)
>>
途径二：%problem2
function dy=problem2(t,y)
dy = zeros(2,1);
dy(1) = 3*y(1)+4*y(2);
dy(2) = -4*y(1)+3*y(2);
[t,y] = ode45('problem2',[0 2],[0 1]);
plot(t,y(:,1),'r',t,y(:,2),'b');
图2
3．假设著名的Lorenz 模型的状态方程表示为
⎪⎩⎪⎨⎧-+-=+-=+-=)
()()()()()()()()()()()(322133223211t x t x t x t x t x t x t x t x t x t x t x t x σρρβ 其中，设28,10,3/8===σρβ。若令其初值为ε===)0(,0)0()0(321x x x ，而ε为机器上可
ε，试用尽可能多的方法来求解该微分方以识别的小常数，如取一个很小的正数10
=
10-
程组。
解：方法一：新建文件lorenzeq.m
源程序： function xdot=lorenzeq(t,x)
xdot=[-8/3*x(1)+x(2)*x(3);...
-10*x(2)+10*x(3);...
-x(1)*x(2)+28*x(2)-x(3)];
再输入命令: t_final=100;
x0=[0;0;1e-10];
[t,x]=ode45('lorenzeq',[0,t_final],x0);
plot(t,x);
figure;
plot3(x(:,1),x(:,2),x(:,3));
axis([10 42 -20 20 -20 25]);
