[分享]三体问题_无聊写的

[huangboat]

function r=sky(t,x)
M=1000;
m1=1;
m2=0.05;
r1=(x(1).^2+x(2).^2+x(3).^2).^1.5;
r2=(x(7).^2+x(8).^2+x(9).^2).^1.5;
r12=((x(1)-x(7)).^2+(x(2)-x(8)).^2+(x(3)-x(9)).^2).^1.5;
r=[x(4);x(5);x(6);-x(1)*M./r1+(x(7)-x(1))*m2./r12;-x(2)*M./r1+(x(8)-x(2))*m2./r12;-x(3)*M./r1+(x(9)-x(3))*m2./r12;...
x(10);x(11);x(12);-x(7)*M./r2+(x(1)-x(7))*m1./r12;-x(8)*M./r2+(x(2)-x(8))*m1./r12;-x(9)*M./r2+(x(3)-x(9))*m1./r12];


clc
clear
w=1e-9*ones(1,12);
q= odeset('RelTol',1e-7,'AbsTol',w);
x0=[2.93;0;0;0;18.45;0;3.1;0;0;0.5;18.70;1];
x1=[2;0;0;0;10;0;10;0;0;0;8;0];
[t,x]=ode45('sky',[0 50],x0,q);
figure(1)
mm=plot3(x(:,1),x(:,2),x(:,3),'r',x(:,7),x(:,8),x(:,9),'b','LineWidth',1)
axis([-5 5 -5 5 -1 1]*2)
xlabel('x')
ylabel('y')
zlabel('z')


figure(2)
r1=(x(:,1).^2+x(:,2).^2+x(:,3).^2).^0.5;
plot(t,r1,'.','MarkerSize',1)
title('m1')


figure(3)
r2=(x(:,7).^2+x(:,8).^2+x(:,9).^2).^0.5;
plot(t,r2,'.','MarkerSize',1)
title('m2')





n=length(t);

mm=plot3(x(:,1),x(:,2),x(:,3),'r.',x(:,7),x(:,8),x(:,9),'b.','LineWidth',1,'MarkerSize',10)
axis([-5 5 -5 5 -1 1])
xlabel('x')
ylabel('y')
zlabel('z')
r1=x(:,1);
r2=x(:,2);
r3=x(:,3);
h1=x(:,7);
h2=x(:,8);
h3=x(:,9);
echo on
set(mm,'erasemode','none'); % none normal xor background
for i=1:n
set(mm,'xdata',r1(i),'ydata',r2(i),'zdata',r3(i),'Color','b','MarkerSize',1)
set(mm,'xdata',h1(i),'ydata',h2(i),'zdata',h3(i),'Color','r','MarkerSize',1)
drawnow;
end
echo off