1条回答 默认 最新
- 技术专家团-Joel 2021-11-13 12:26关注
你好,可以参考为你写的代码(如有帮助还望题主给个宝贵的采纳支持一下答主哟):
微分方程:Matlab代码(有动图)
% x(1)对应r,x(2)对应dr/dt,x(3)对应theta,x(4)对应dtheta/dt M=5.965e24; %地球质量 G=6.67e-11; %万有引力常数 odefun = @(t,x)[ x(2); x(1)*x(4)^2-G*M/x(1)^2; x(4); -2*x(2)*x(4)/x(1); ]; r0 = 384403.9e3; %地月距离 dr0 = 0;%初始径速度 theta0 = 0;%初始角度 dtheta0 = 1020/r0;%初始角速度 x0 = [r0,dr0,theta0,dtheta0]; [t,y] = ode45(odefun,[0,3600*24*30*3],x0);%模拟三个月 r = y(:,1); theta = y(:,3); xp = r.*cos(theta); yp = r.*sin(theta); % 画图 for i = 1:1:numel(xp) figure(1);clf plot(xp(1:i)/1e7,yp(1:i)/1e7,'b-') hold on; plot(0,0,'bo','markersize',15,'markerfacecolor','b') plot(xp(i)/1e7,yp(i)/1e7,'mo','markersize',5,'markerfacecolor','m') axis([-1,1,-1,1]*1.2*r0/1e7) xlabel('X(万公里)');ylabel('Y(万公里)') text(xp(i)/1e7+3,yp(i)/1e7+3,'月球') text(3,3,'地球') title(['第',num2str(t(i)/3600/24),'天']) axis equal if(mod(i,3600)==1) pause(0.00001) end end
图片效果:
本回答被题主选为最佳回答 , 对您是否有帮助呢?解决 2无用