Three body problem simulation for MatLab.
% Figure8
p1 = 0.347111;
p2 = 0.532728;
% Butterfly
pp1 = 0.306893;
pp2 = 0.125507;
% Butterfly 4
ppppp1 = 0.350112;
ppppp2 = 0.079339;
% MothIII
ppp1 = 0.383444;
ppp2 = 0.377364;
% Ying Yang
pppp1 = 0.416822;
pppp2 = 0.330333;
% Dragonfly
pppppp1 = 0.080584;
pppppp2 = 0.588836;
% Yarn
ppppppp1 = 0.55906;
ppppppp2 = 0.34919;
% Goggles
pppppppp1 = 0.083300;
pppppppp2 = 0.127889;
m1 = 1;
m2 = 1;
m3 = 1;
x1_0 = -1;
y1_0 = 0;
z1_0 = 0;
x2_0 = 1;
y2_0 = 0;
z2_0 = 0;
x3_0 = 0;
y3_0 = 0;
z3_0 = 0;
vx1_0 = pp1;
vy1_0 = pp2;
vz1_0 = 0.1;
vx2_0 = pp1;
vy2_0 = pp2;
vz2_0 = -0.1;
vx3_0 = -2pp1;
vy3_0 = -2pp2;
vz3_0 = 0;
t0 = 0;
tn = 50;
r = @ (a,b,c) sqrt(a^2 + b^2 + c^2);
% u = [x,y,z,vx,vy,vz]
u01 = [x1_0, y1_0, z1_0];
u02 = [x2_0, y2_0, z2_0];
u03 = [x3_0, y3_0, z3_0];
u04 = [vx1_0, vy1_0, vz1_0];
u05 = [vx2_0, vy2_0, vz2_0];
u06 = [vx3_0, vy3_0, vz3_0];
u0 = [u01, u02, u03, u04, u05, u06];
M = @(t,u) [u(10); %First Body X 1
u(11); %First Body Y
u(12); %First Body Z
u(13); %Second Body X 4
u(14); %Second Body Y
u(15); %Second Body Z
u(16); %Third Body X 7
u(17); %Third Body Y
u(18); %Third Body Z
-m2(u(1)-u(4))/(r(u(1)-u(4),u(2)-u(5),u(3)-u(6)))^3 – m3(u(1)-u(7))/(r(u(1)-u(7),u(2)-u(8),u(3)-u(9)))^3;
-m2(u(2)-u(5))/(r(u(1)-u(4),u(2)-u(5),u(3)-u(6)))^3 – m3(u(2)-u(8))/(r(u(1)-u(7),u(2)-u(8),u(3)-u(9)))^3;
-m2(u(3)-u(6))/(r(u(1)-u(4),u(2)-u(5),u(3)-u(6)))^3 – m3(u(3)-u(9))/(r(u(1)-u(7),u(2)-u(8),u(3)-u(9)))^3;
-m3(u(4)-u(7))/(r(u(4)-u(7),u(5)-u(8),u(6)-u(9)))^3 – m1(u(4)-u(1))/(r(u(4)-u(1),u(5)-u(2),u(6)-u(3)))^3;
-m3(u(5)-u(8))/(r(u(4)-u(7),u(5)-u(8),u(6)-u(9)))^3 – m1(u(5)-u(2))/(r(u(4)-u(1),u(5)-u(2),u(6)-u(3)))^3;
-m3(u(6)-u(9))/(r(u(4)-u(7),u(5)-u(8),u(6)-u(9)))^3 – m1(u(6)-u(3))/(r(u(4)-u(1),u(5)-u(2),u(6)-u(3)))^3;
-m1(u(7)-u(1))/(r(u(7)-u(1),u(8)-u(2),u(9)-u(3)))^3 – m2(u(7)-u(4))/(r(u(7)-u(4),u(8)-u(5),u(9)-u(6)))^3;
-m1(u(8)-u(2))/(r(u(7)-u(1),u(8)-u(2),u(9)-u(3)))^3 – m2(u(8)-u(5))/(r(u(7)-u(4),u(8)-u(5),u(9)-u(6)))^3;
-m1(u(9)-u(3))/(r(u(7)-u(1),u(8)-u(2),u(9)-u(3)))^3 – m2(u(9)-u(6))/(r(u(7)-u(4),u(8)-u(5),u(9)-u(6)))^3;];
options = odeset(‘RelTol’,1e-6,’AbsTol’,1e-8);
sol = ode45(M,t0:tn,u0,options);
% sol2 = ode45(M,t0:tn,u0)
% fplot(@(x)deval(sol,x,1), [t0, tn]) % x
% xlabel(‘t’)
% ylabel(‘x’)
% fplot(@(x)deval(sol,x,2), [t0, tn]) % y
% xlabel(‘t’)
% ylabel(‘y’)
% fplot(@(x)deval(sol,x,3), [t0, tn]) % vx
% xlabel(‘t’)
% ylabel(‘vx’)
% fplot(@(x)deval(sol,x,4), [t0, tn]) % vy
% xlabel(‘t’)
% ylabel(‘vy’)
x = linspace(t0,tn,5000);
body1x = deval(sol,x,1); %1x
body1y = deval(sol,x,2); %1y
body1z = deval(sol,x,3); %1z
body2x = deval(sol,x,4); %2x
body2y = deval(sol,x,5); %2y
body2z = deval(sol,x,6); %2z
body3x = deval(sol,x,7); %3x
body3y = deval(sol,x,8); %3y
body3z = deval(sol,x,9); %3z
% body1vx = deval(sol,x,10); %1vx
% body1vy = deval(sol,x,11); %1vy
% body1vz = deval(sol,x,12); %1vz
% body2vx = deval(sol,x,13); %2vx
% body2vy = deval(sol,x,14); %2vy
% body2vz = deval(sol,x,15); %2vz
% body3vx = deval(sol,x,16); %3vx
% body3vy = deval(sol,x,17); %3vy
% body3vz = deval(sol,x,18); %3vz
Lx = max([body1x,body2x,body3x]);
Ly = max([body1y,body2y,body3y]);
Lz = max([body1z,body2z,body3z]);
lx = min([body1x,body2x,body3x]);
ly = min([body1y,body2y,body3y]);
lz = min([body1z,body2z,body3z]);
% figure(‘Name’,’Body 1′);
% plot3(body1x,body1y,body1z,’color’,[0 0.4470 0.7410]); %x vs y vs z
% xlabel(‘x’);
% ylabel(‘y’);
% zlabel(‘z’);
% view(45,25);
% grid on;
% legend(strcat(‘m1 = ‘, num2str(m1)),’location’,’best’);
%
% figure(‘Name’,’Body 2′);
% plot3(body2x,body2y,body2z,’color’,[0.8500 0.3250 0.0980]); %vx vs vy vs vz
% xlabel(‘x’);
% ylabel(‘y’);
% zlabel(‘z’);
% view(45,25);
% grid on;
% legend(strcat(‘m2 = ‘, num2str(m2)),’location’,’best’);
%
% figure(‘Name’,’Body 3′);
% plot3(body3x,body3y,body3z,’color’,[0.9290 0.6940 0.1250]); %vx vs vy vs vz
% xlabel(‘x’);
% ylabel(‘y’);
% zlabel(‘z’);
% view(45,25);
% grid on;
% legend(strcat(‘m3 = ‘, num2str(m3)),’location’,’best’);
%
% figure(‘Name’,’Body Paths’);
% plot3(body1x,body1y,body1z,body2x,body2y,body2z,body3x,body3y,body3z);
% view(45,45);
% grid on;
% set(gca,’XLim’,[lx-0.02 Lx+0.02],’YLim’,[ly-0.02 Ly+0.02],’ZLim’,[lz-0.02 Lz+0.02]);
figure(‘Name’,’3 Body Problem’);
curve1 = animatedline(‘color’,[0 0.4470 0.7410]);
hold on
curve2 = animatedline(‘color’,[0.8500 0.3250 0.0980]);
curve3 = animatedline(‘color’,[0.9290 0.6940 0.1250]);
curve4 = animatedline(‘Marker’,’.’,’MarkerSize’,10m1+1,’MaximumNumPoints’,1,’color’,[0 0.4470 0.7410]);
curve5 = animatedline(‘Marker’,’.’,’MarkerSize’,10m2+1,’MaximumNumPoints’,1,’color’,[0.8500 0.3250 0.0980]);
curve6 = animatedline(‘Marker’,’.’,’MarkerSize’,10m3+1,’MaximumNumPoints’,1,’color’,[0.9290 0.6940 0.1250]);
set(gca,’XLim’,[lx-0.02 Lx+0.02],’YLim’,[ly-0.02 Ly+0.02],’ZLim’,[lz-0.02 Lz+0.02]);
view(45,25);
grid on;
xlabel(‘x’);
ylabel(‘y’);
zlabel(‘z’);
rotate3d;
for i = 1:length(x)
%view(90i/length(x),90*i/length(x));
addpoints(curve1,body1x(i),body1y(i),body1z(i));
addpoints(curve2,body2x(i),body2y(i),body2z(i));
addpoints(curve3,body3x(i),body3y(i),body3z(i));
addpoints(curve4,body1x(i),body1y(i),body1z(i));
addpoints(curve5,body2x(i),body2y(i),body2z(i));
addpoints(curve6,body3x(i),body3y(i),body3z(i));
drawnow limitrate
pause(0.0001)
end
legend(strcat(‘m1 = ‘, num2str(m1)),strcat(‘m2 = ‘, num2str(m2)),strcat(‘m3 = ‘, num2str(m3)),’location’,’best’);
