# Parker-Sochacki 3-Body Gravitation Problem

Greetings! In previous blogs I applied the Parker-Sochack method to two versions of the 2-Body Gravitation Problem. One involved two moving masses where I showed the benefits of a higher order method. This was done in 2D and 3D. The other involved a moving mass and stationary point where I showed interesting geometry patterns that could be used as particle effects for a game or film.

In this blog, I expand the problem to a 3-Body Gravitation problem in 2D solved with the Parker-Sochacki method. I will show how I implemented the math, code in Matlab, output of different initial conditions, and discuss benefits of using a higher order method. Now that I expanded the 2-Body problem to a 3-Body problem, this should also be a useful setup to generalize to an n-body problem where n is greater than or equal to 3.

I would like to make two points before showing the implementation. The first point is that my method is by no means the optimal solution. I treated the the separation distance between mass 1 and mass 2 differently than the separation distance between mass 2 and mass 1 when I could have used the same variable and been careful with signs. I did this to make my code more readable and make a point a little later in this blog. My second point is this. Higher order methods are well suited for close encounters. I mean that a higher order method does well when objects are close together. In this case, when objects were close together under gravitational forces, I needed to increase the order of Parker-Sochacki at a time step of 0.001 otherwise the objects would “blow up” diverging from each other. Not even the accuracy of Runge-Kutta 4th order could handle this simulation at a time step of 0.001 with the initial conditions I played with. This means that Parker-Sochacki is well suited for any n-body type simulation when the points are close together.

The first thing I will show are the First Order Differential Equations that make up the n-body problem under gravity in 3-Dimensions:

We have particle “j” with x,y, and z components. The derivative of position for particle j is particle j’s velocity for x, y, and z. The derivative of velocity is a summation where we sum over all of the neighboring particles(This blog will not cover simplifications to this). We have a different particle “k” THAT IS NOT EQUAL TO “j”(look at the math and see that we would have zero in the denominator creating problems) and we take the difference of position “k” minus position “j” multiplied by the mass of particle “k” and divide it by the norm of “j” minus “k” (note the order switches). Then, while “j” is the same particle, we now change “k” to another particle and continue this for all particles. This is all multiplied by the universal gravitation constant G.

Now I will do this for 3 particles in 2-D cutting out the z-component. I will show the ODEs for all three particles calling them particle 1, particle 2, and particle 3. I’m expanding the summation sign and switching signs(now you see minus signs in front).

To apply the Parker-Sochacki method, we introduce auxiliary variables in the denominator. I will introduce 6 auxiliary variables:

Note that auxiliary variable x and auxiliary variable y are not x and y positions. I’m running out of variable choices when I am done with this. With these auxiliary variables, our ODEs for Velocity become:

The next step is to take derivatives for s, u, v, w, x, and y.

This simplifies to:

When multiplying terms, we won’t do more than 2 at once using cauchy products. This means that we need to introduce another 6 auxiliary variables:

This gets plugged back into s’, u’, v’, w’, x’, and y’:

We now have everything we need to create a system of First Order ODEs

The goal with this system of ODEs is to have a Maclaurin series for the position, velocity, and separation variables. I will show what the recurrence relationship looks like. You will see how we handle any cubic terms s,u,v,w,x,y and see that we do not need derivatives for a,b,c,d,e, and f. When looking at the number subscripts for position and velocity, note that I am referring to “particle 1”, “particle 2”, and “particle 3.” That does not go towards calculating the i+1 term. Given an initial set of positions and velocities, we can create an initial set of conditions for all of our variables. We then use those values to create Maclaurin coefficients up to any order desired:

I’m copying and pasting from writing the formulas into Latex. Note that I created s2,u2,v2,w2,x2,y2,s3,u3,v3,w3,x3, and y3 to handle the cube terms. Note that all of those variables plus a,b,c,d,e, and f are computed at coefficient i instead of i+1. I do not need their derivatives as states earlier.

Given a time step of your choice, t, we create a Maclaurin series for the following variables:

I am not creating a Maclaurin series for any term that computes i instead of i+1 because those initial conditions are created from these variables anyway. Once this Maclaurin series is calculated, these variables are updated becoming the new initial conditions for the next iteration.

Next is the code. I wrote this in Matlab.

% June 6, 2016
% This program uses the Parker-Sochacki Method to solve the 3-Body probolem
% in 2 dimensions

% clear everything
close all;
clear all;
clc;

% Parameters
G = 1;
m1 = 0.6;
m2 = 0.1001;
m3 = 0.00101;

% Declare variables for Parker-Sochacki Implementation
order = 30; % order for power series
t = 0.001; % time step
frames = 2125; % number of iterations

% Create arrays for all Maclaurin coefficients
x_pos1 = zeros(1,order-1);
y_pos1 = zeros(1,order-1);
vx1 = zeros(1,order-1);
vy1 = zeros(1,order-1);
x_pos2 = zeros(1,order-1);
y_pos2 = zeros(1,order-1);
vx2 = zeros(1,order-1);
vy2 = zeros(1,order-1);
x_pos3 = zeros(1,order-1);
y_pos3 = zeros(1,order-1);
vx3 = zeros(1,order-1);
vy3 = zeros(1,order-1);
s = zeros(1,order-1);
s2 = zeros(1,order-1);
s3 = zeros(1,order-1);
a = zeros(1,order-1);
u = zeros(1,order-1);
u2 = zeros(1,order-1);
u3 = zeros(1,order-1);
b = zeros(1,order-1);
v = zeros(1,order-1);
v2 = zeros(1,order-1);
v3 = zeros(1,order-1);
c = zeros(1,order-1);
w = zeros(1,order-1);
w2 = zeros(1,order-1);
w3 = zeros(1,order-1);
d = zeros(1,order-1);
x = zeros(1,order-1);
x2 = zeros(1,order-1);
x3 = zeros(1,order-1);
e = zeros(1,order-1);
y = zeros(1,order-1);
y2 = zeros(1,order-1);
y3 = zeros(1,order-1);
f = zeros(1,order-1);

% initial mass 1 conditions
x_pos1(1) = 0;
vx1(1) = 0.60;
y_pos1(1) = 1.95;
vy1(1) = 0.5;

% initial mass 2 conditions
x_pos2(1) = 0;
vx2(1) = 2.828;
y_pos2(1) = 2;
vy2(1) = 0;

% initial mass 3 conditions
x_pos3(1) = 0;
vx3(1) = 3.4;
y_pos3(1) = 2.1;
vy3(1) = 0;

% initial auxiliary variable 1 conditions 1-2
s(1) = 1./sqrt((x_pos1(1) – x_pos2(1)).^2 + (y_pos1(1) – y_pos2(1)).^2);
s2(1) = s(1)*s(1);
s3(1) = s2(1)*s(1);
a(1) = (x_pos1(1) – x_pos2(1))*(vx1(1) – vx2(1)) + (y_pos1(1) – y_pos2(1))*(vy1(1)-vy2(1));

% initial auxiliary variable 2 conditions 2-1
u(1) = 1./sqrt((x_pos2(1) – x_pos1(1)).^2 + (y_pos2(1) – y_pos1(1)).^2);
u2(1) = u(1)*u(1);
u3(1) = u2(1)*u(1);
b(1) = (x_pos2(1) – x_pos1(1))*(vx2(1) – vx1(1)) + (y_pos2(1) – y_pos1(1))*(vy2(1)-vy1(1));

% initial auxiliary variable 3 conditions 1-3
v(1) = 1./sqrt((x_pos1(1) – x_pos3(1)).^2 + (y_pos1(1) – y_pos3(1)).^2);
v2(1) = v(1)*v(1);
v3(1) = v2(1)*v(1);
c(1) = (x_pos1(1) – x_pos3(1))*(vx1(1) – vx3(1)) + (y_pos1(1) – y_pos3(1))*(vy1(1)-vy3(1));

% initial auxiliary variable 4 conditions 3-1
w(1) = 1./sqrt((x_pos3(1) – x_pos1(1)).^2 + (y_pos3(1) – y_pos1(1)).^2);
w2(1) = w(1)*w(1);
w3(1) = w2(1)*w(1);
d(1) = (x_pos3(1) – x_pos1(1))*(vx3(1) – vx1(1)) + (y_pos3(1) – y_pos1(1))*(vy3(1)-vy1(1));

% initial auxiliary variable 5 conditions 2-3
x(1) = 1./sqrt((x_pos2(1) – x_pos3(1)).^2 + (y_pos2(1) – y_pos3(1)).^2);
x2(1) = x(1)*x(1);
x3(1) = x2(1)*x(1);
e(1) = (x_pos2(1) – x_pos3(1))*(vx2(1) – vx3(1)) + (y_pos2(1) – y_pos3(1))*(vy2(1)-vy3(1));

% initial auxiliary variable 6 conditions 3-2
y(1) = 1./sqrt((x_pos3(1) – x_pos2(1)).^2 + (y_pos3(1) – y_pos2(1)).^2);
y2(1) = y(1)*y(1);
y3(1) = y2(1)*y(1);
f(1) = (x_pos3(1) – x_pos2(1))*(vx3(1) – vx2(1)) + (y_pos3(1) – y_pos2(1))*(vy3(1)-vy2(1));

% mass 1 power series sum
x_pos1_sum = zeros(1,frames);
y_pos1_sum = zeros(1,frames);
vx1_sum = zeros(1,frames);
vy1_sum = zeros(1, frames);

% mass 2 power series sum
x_pos2_sum = zeros(1,frames);
y_pos2_sum = zeros(1,frames);
vx2_sum = zeros(1,frames);
vy2_sum = zeros(1, frames);

% mass 3 power series sum
x_pos3_sum = zeros(1,frames);
y_pos3_sum = zeros(1,frames);
vx3_sum = zeros(1,frames);
vy3_sum = zeros(1, frames);

% Auxiliary Variable sums
s_sum = zeros(1, frames);
s2_sum = zeros(1, frames);
s3_sum = zeros(1, frames);
a_sum = zeros(1, frames);
u_sum = zeros(1, frames);
u2_sum = zeros(1, frames);
u3_sum = zeros(1, frames);
b_sum = zeros(1, frames);
v_sum = zeros(1, frames);
v2_sum = zeros(1, frames);
v3_sum = zeros(1, frames);
c_sum = zeros(1, frames);
w_sum = zeros(1, frames);
w2_sum = zeros(1, frames);
w3_sum = zeros(1, frames);
d_sum = zeros(1, frames);
x_sum = zeros(1, frames);
x2_sum = zeros(1, frames);
x3_sum = zeros(1, frames);
e_sum = zeros(1, frames);
y_sum = zeros(1, frames);
y2_sum = zeros(1, frames);
y3_sum = zeros(1, frames);
f_sum = zeros(1, frames);

% Give initial value sums
x_pos1_sum(1) = x_pos1(1);
y_pos1_sum(1) = y_pos1(1);
vx1_sum(1) = vx1(1);
vy1_sum(1) = vy1(1);
x_pos2_sum(1) = x_pos2(1);
y_pos2_sum(1) = y_pos2(1);
vx2_sum(1) = vx2(1);
vy2_sum(1) = vy2(1);
x_pos3_sum(1) = x_pos3(1);
y_pos3_sum(1) = y_pos3(1);
vx3_sum(1) = vx3(1);
vy3_sum(1) = vy3(1);
s_sum(1) = s(1);
u_sum(1) = u(1);
v_sum(1) = v(1);
w_sum(1) = w(1);
x_sum(1) = x(1);
y_sum(1) = y(1);

% This for loop updates the frames
for j= 2:frames

% Calculate Coefficients
for iterate = 1 : order-1

% Setup Cauchy Product
% Velocity
cp_vx1 = 0;
cp_vy1 = 0;
cp_vx2 = 0;
cp_vy2 = 0;
cp_vx3 = 0;
cp_vy3 = 0;

% Auxiliary Variables Cauchy Product 1-2
cp_s = 0;
cp_s2 = 0;
cp_s3 = 0;
cp_a = 0;

% Auxiliary Variables Cauchy Product 2-1
cp_u = 0;
cp_u2 = 0;
cp_u3 = 0;
cp_b = 0;

% Auxiliary Variables Cauchy Product 1-3
cp_v = 0;
cp_v2 = 0;
cp_v3 = 0;
cp_c = 0;

% Auxiliary Variables Cauchy Product 3-1
cp_w = 0;
cp_w2 = 0;
cp_w3 = 0;
cp_d = 0;

% Auxiliary Variables Cauchy Product 2-3
cp_x = 0;
cp_x2 = 0;
cp_x3 = 0;
cp_e = 0;

% Auxiliary Variables Cauchy Product 3-2
cp_y = 0;
cp_y2 = 0;
cp_y3 = 0;
cp_f = 0;

% Cauchy Products s2, u2, v2, w2, x2, y2
for k = 1:iterate

cp_s2 = cp_s2 + (s(k).*s(iterate+1-k));
cp_u2 = cp_u2 + (u(k).*u(iterate+1-k));
cp_v2 = cp_v2 + (v(k).*v(iterate+1-k));
cp_w2 = cp_w2 + (w(k).*w(iterate+1-k));
cp_x2 = cp_x2 + (x(k).*x(iterate+1-k));
cp_y2 = cp_y2 + (y(k).*y(iterate+1-k));

end

% Cauchy Product update s2, u2, v2, w2, x2, y2
s2(iterate) = cp_s2; % s2
u2(iterate) = cp_u2; % u2
v2(iterate) = cp_v2; % v2
w2(iterate) = cp_w2; % w2
x2(iterate) = cp_x2; % x2
y2(iterate) = cp_y2; % y2

% Cauchy Products s3
for k = 1:iterate

cp_s3 = cp_s3 + (s2(k).*s(iterate+1-k));
cp_u3 = cp_u3 + (u2(k).*u(iterate+1-k));
cp_v3 = cp_v3 + (v2(k).*v(iterate+1-k));
cp_w3 = cp_w3 + (w2(k).*w(iterate+1-k));
cp_x3 = cp_x3 + (x2(k).*x(iterate+1-k));
cp_y3 = cp_y3 + (y2(k).*y(iterate+1-k));

end

% Cauchy Product update s3, u3, v3, w3, x3, y3
s3(iterate) = cp_s3; % s3
u3(iterate) = cp_u3; % u3
v3(iterate) = cp_v3; % v3
w3(iterate) = cp_w3; % w3
x3(iterate) = cp_x3; % x3
y3(iterate) = cp_y3; % y3

% Cauchy Products a, b, c, d, e, f
for k = 1:iterate

cp_a = cp_a + ((x_pos1(k) – x_pos2(k)).*(vx1((iterate+1)-k)-vx2((iterate+1)-k)) + (y_pos1(k) – y_pos2(k)).*((vy1((iterate+1)-k)-vy2((iterate+1)-k))));
cp_b = cp_b + ((x_pos2(k) – x_pos1(k)).*(vx2((iterate+1)-k)-vx1((iterate+1)-k)) + (y_pos2(k) – y_pos1(k)).*((vy2((iterate+1)-k)-vy1((iterate+1)-k))));
cp_c = cp_c + ((x_pos1(k) – x_pos3(k)).*(vx1((iterate+1)-k)-vx3((iterate+1)-k)) + (y_pos1(k) – y_pos3(k)).*((vy1((iterate+1)-k)-vy3((iterate+1)-k))));
cp_d = cp_d + ((x_pos3(k) – x_pos1(k)).*(vx3((iterate+1)-k)-vx1((iterate+1)-k)) + (y_pos3(k) – y_pos1(k)).*((vy3((iterate+1)-k)-vy1((iterate+1)-k))));
cp_e = cp_e + ((x_pos2(k) – x_pos3(k)).*(vx2((iterate+1)-k)-vx3((iterate+1)-k)) + (y_pos2(k) – y_pos3(k)).*((vy2((iterate+1)-k)-vy3((iterate+1)-k))));
cp_f = cp_f + ((x_pos3(k) – x_pos2(k)).*(vx3((iterate+1)-k)-vx2((iterate+1)-k)) + (y_pos3(k) – y_pos2(k)).*((vy3((iterate+1)-k)-vy2((iterate+1)-k))));

end

% Cauchy Product update a, b, c, d, e, f
a(iterate) = cp_a; % a
b(iterate) = cp_b; % b
c(iterate) = cp_c; % c
d(iterate) = cp_d; % d
e(iterate) = cp_e; % e
f(iterate) = cp_f; % f

% Cauchy Products for separation auxiliary variable s, u, v, w, x, y
for k =1:iterate

cp_s = cp_s + (s3(k)*a((iterate+1)-k));
cp_u = cp_u + (u3(k)*b((iterate+1)-k));
cp_v = cp_v + (v3(k)*c((iterate+1)-k));
cp_w = cp_w + (w3(k)*d((iterate+1)-k));
cp_x = cp_x + (x3(k)*e((iterate+1)-k));
cp_y = cp_y + (y3(k)*f((iterate+1)-k));

end

% Auxiliary Variable separation update
s(iterate+1) = -cp_s/iterate; % s
u(iterate+1) = -cp_u/iterate; % u
v(iterate+1) = -cp_v/iterate; % v
w(iterate+1) = -cp_w/iterate; % w
x(iterate+1) = -cp_x/iterate; % x
y(iterate+1) = -cp_y/iterate; % y

% Position for mass 1
x_pos1(iterate+1) = vx1(iterate)/iterate;
y_pos1(iterate+1) = vy1(iterate)/iterate;

% Position for mass 2
x_pos2(iterate+1) = vx2(iterate)/iterate;
y_pos2(iterate+1) = vy2(iterate)/iterate;

% Position for mass 3
x_pos3(iterate+1) = vx3(iterate)/iterate;
y_pos3(iterate+1) = vy3(iterate)/iterate;

% Cauchy Products for Velocity
for k =1:iterate

% Velocity x for mass 1
cp_vx1 = cp_vx1 – (G*m2*(x_pos1(k)-x_pos2(k))*s3((iterate+1)-k));
cp_vx1 = cp_vx1 – (G*m3*(x_pos1(k)-x_pos3(k))*v3((iterate+1)-k));

% Velocity y for mass 1
cp_vy1 = cp_vy1 – (G*m2*(y_pos1(k)-y_pos2(k))*s3((iterate+1)-k));
cp_vy1 = cp_vy1 – (G*m3*(y_pos1(k)-y_pos3(k))*v3((iterate+1)-k));

% Velocity x for mass 2
cp_vx2 = cp_vx2 – (G*m1*(x_pos2(k)-x_pos1(k))*u3((iterate+1)-k));
cp_vx2 = cp_vx2 – (G*m3*(x_pos2(k)-x_pos3(k))*x3((iterate+1)-k));

% Velocity y for mass 2
cp_vy2 = cp_vy2 – (G*m1*(y_pos2(k)-y_pos1(k))*u3((iterate+1)-k));
cp_vy2 = cp_vy2 – (G*m3*(y_pos2(k)-y_pos3(k))*x3((iterate+1)-k));

% Velocity x for mass 3
cp_vx3 = cp_vx3 – (G*m1*(x_pos3(k)-x_pos1(k))*w3((iterate+1)-k));
cp_vx3 = cp_vx3 – (G*m2*(x_pos3(k)-x_pos2(k))*y3((iterate+1)-k));

% Velocity y for mass 3
cp_vy3 = cp_vy3 – (G*m1*(y_pos3(k)-y_pos1(k))*w3((iterate+1)-k));
cp_vy3 = cp_vy3 – (G*m2*(y_pos3(k)-y_pos2(k))*y3((iterate+1)-k));

end

% Cauchy Product update vx1, vy1, vx2, vy2, vx3, vy3
% Velocity for mass 1
vx1(iterate+1) = cp_vx1/iterate;
vy1(iterate+1) = cp_vy1/iterate;

% Velocity for mass 2
vx2(iterate+1) = cp_vx2/iterate;
vy2(iterate+1) = cp_vy2/iterate;

% Velocity for mass 3
vx3(iterate+1) = cp_vx3/iterate;
vy3(iterate+1) = cp_vy3/iterate;

end

%
% % Setup sums
x1_temp_sum = 0;
y1_temp_sum = 0;
vx1_temp_sum = 0;
vy1_temp_sum = 0;
x2_temp_sum = 0;
y2_temp_sum = 0;
vx2_temp_sum = 0;
vy2_temp_sum = 0;
x3_temp_sum = 0;
y3_temp_sum = 0;
vx3_temp_sum = 0;
vy3_temp_sum = 0;
s_temp_sum = 0;
u_temp_sum = 0;
v_temp_sum = 0;
w_temp_sum = 0;
x_temp_sum = 0;
y_temp_sum = 0;

% Compute Power series
for i = 1:order-1

x1_temp_sum = x1_temp_sum + (x_pos1(i)*t^(i-1));
y1_temp_sum = y1_temp_sum + (y_pos1(i)*t^(i-1));
vx1_temp_sum = vx1_temp_sum + (vx1(i)*t^(i-1));
vy1_temp_sum = vy1_temp_sum + (vy1(i)*t^(i-1));
x2_temp_sum = x2_temp_sum + (x_pos2(i)*t^(i-1));
y2_temp_sum = y2_temp_sum + (y_pos2(i)*t^(i-1));
vx2_temp_sum = vx2_temp_sum + (vx2(i)*t^(i-1));
vy2_temp_sum = vy2_temp_sum + (vy2(i)*t^(i-1));
x3_temp_sum = x3_temp_sum + (x_pos3(i)*t^(i-1));
y3_temp_sum = y3_temp_sum + (y_pos3(i)*t^(i-1));
vx3_temp_sum = vx3_temp_sum + (vx3(i)*t^(i-1));
vy3_temp_sum = vy3_temp_sum + (vy3(i)*t^(i-1));
s_temp_sum = s_temp_sum + (s(i)*t^(i-1));
u_temp_sum = u_temp_sum + (u(i)*t^(i-1));
v_temp_sum = v_temp_sum + (v(i)*t^(i-1));
w_temp_sum = w_temp_sum + (w(i)*t^(i-1));
x_temp_sum = x_temp_sum + (x(i)*t^(i-1));
y_temp_sum = y_temp_sum + (y(i)*t^(i-1));

end

% % Put in power series to store in array
x_pos1_sum(j) = x1_temp_sum;
y_pos1_sum(j) = y1_temp_sum;
vx1_sum(j) = vx1_temp_sum;
vy1_sum(j) = vy1_temp_sum;
x_pos2_sum(j) = x2_temp_sum;
y_pos2_sum(j) = y2_temp_sum;
vx2_sum(j) = vx2_temp_sum;
vy2_sum(j) = vy2_temp_sum;
x_pos3_sum(j) = x3_temp_sum;
y_pos3_sum(j) = y3_temp_sum;
vx3_sum(j) = vx3_temp_sum;
vy3_sum(j) = vy3_temp_sum;
s_sum(j) = s_temp_sum;
u_sum(j) = u_temp_sum;
v_sum(j) = v_temp_sum;
w_sum(j) = w_temp_sum;
x_sum(j) = x_temp_sum;
y_sum(j) = y_temp_sum;

%
%
% % reset initial conditions
x_pos1(1) = x1_temp_sum;
y_pos1(1) = y1_temp_sum;
vx1(1) = vx1_temp_sum;
vy1(1) = vy1_temp_sum;
x_pos2(1) = x2_temp_sum;
y_pos2(1) = y2_temp_sum;
vx2(1) = vx2_temp_sum;
vy2(1) = vy2_temp_sum;
x_pos3(1) = x3_temp_sum;
y_pos3(1) = y3_temp_sum;
vx3(1) = vx3_temp_sum;
vy3(1) = vy3_temp_sum;
s(1) = s_temp_sum;
u(1) = u_temp_sum;
v(1) = v_temp_sum;
w(1) = w_temp_sum;
x(1) = x_temp_sum;
y(1) = y_temp_sum;

end

% Plot Results
plot(x_pos1_sum,y_pos1_sum,x_pos2_sum,y_pos2_sum, x_pos3_sum, y_pos3_sum);
legend(‘mass 1’, ‘mass 2’, ‘mass 3’);
title(‘Parker-Sochacki 3-Body Implementation’);

Given different masses, initial positions, and initial velocities, we got the following pattern:

In blue is the heaviest mass followed by the red and then orange.

This implementation went for 2125 frames at a time step of 0.001 and order 30(for Matlab meaning this is order 29 in other languages like C++). Look what happens when I lower this method to the accuracy of Runge-Kutta 4 at the same time step:

The second mass is not able to handle “close encounters” properly. I stated in the beginning that higher order methods do well when the objects are close together.

This simulation, given its time step and number of frames, only goes for 2 seconds. If I keep this to RK4 accuracy, same time step, but change it to 3000 frames for 3 seconds, you see that mass 1(blue) and mass 2(red) are not working properly:

Mass 2 moves in a straight line a little after 2 seconds. RK4 would not be able to handle this simulation unless you lower the time step and this excludes methods like Euler or Verlet at this time step. They blew up right away.

There is another note about order. Generally, a higher order leads to a more accurate solution but not all of the time. For this simulation, 4th order was the first setup to have a pattern that did not blow up. When I increased the implementation to order 5 and then 6 the simulation looked more stable. However, when I increased the implementation to orders 7 and 8, the simulation blew up. Then order 9 looked the better than the previous orders. It was not until I got to order 11 that I had a simulation that looked like the first image you see. For this implementation, at this time step, order 11 was the smallest order to achieve a stable and accurate result.

I want to show another example where one mass is significantly larger than the other 2. This sort of simulates the sun, moon, and earth. For this, I increased the order to 100, time step is the same at 0.001, there are 10000 frames, mass 1 is 16, mass 2 is 0.1001, and mass 3 is 0.00101.:

If you look closely, you see that mass 1(blue) barely moves but it moves. Then mass 2(red-orange) revolves around mass 1. Mass 3(yellow-orange) moves around mass 2.

This is a simulation that can handle a lower order. All conditions the same, when I changed this simulation to the accuracy of RK4, the result looked about the same. However, I could not increase the time step unless I increased the order.

CLOSING THOUGHTS: There are several ways to optimize this simulation. The distance between particle 1 and particle 2 and the distance between particle 2 and particle 1 are the same. I did not do that. This means that I could have cut separation auxiliary variables in half. We could have done with only s,u,v,s2,u2,v2,a,b, and c. This would speed up the simulation with the same result.Also, there is a formula to take the power of a power series shown below:

Doing this would cut s2,u2,v2,w2,x2, and y2 variables cutting cauchy products. This would change how s3,u3,v3,w3,x3,and y3 are calculated. I do not know if this would cut time overall. Then there is the consideration on how to update the auxiliary variables. Pending on the order of the Maclaurin series, is it better to update s,u,v,w,x,y with a Maclaurin series or calculating them as norms like we do in the very beginning?

This post is a setup to apply Parker-Sochacki to the n-body problem. Again, from the different initial conditions I played with, Parker-Sochacki should be used when objects are “close enough together.”

Thank you for taking the time to go through this!

Sincerely,