There are previous blog posts I have written about applying the Parker-Sochacki method to the 2 and 3-body gravity problems. I discussed results(showing benefits of a higher order method and how Parker-Sochacki does well against close encounters(bodies that are very close).

This blog picks up from the previous blogs. I still take the power of a power series so I skip auxiliary variable s2.

At long last, I have generalized the 3-body problem into an N-body problem. The code posted below can be modified to as many bodies as you want. All you need to do is change num_of_bodies to the number of bodies in your simulation, add, remove, or change the initial conditions, and then adjust how you plot the final solution.

First, I will show the Ordinary Differential Equations again from a previous post, then I will show the code, then I will discuss what I did involving matrices, then I will show output for number of bodies set to 2, 3, and 5.

In 3D, these are the equations of motion to update the position and velocity of particle j:

I implemented the n-body problem in 2D. There is a previous blog where I implemented a 3D version of the 2-body problem in Houdini.

For the 3-body problem, I had an auxiliary variable represent particle 1 to particle 2, particle 1 to particle 3, and particle 2 to particle 3. It was tedious being hard coded. I used matrices in Matlab to generalize this.

I use the same notation that is done in the published Parker-Sochacki n-body papers out of respect. They were the geniuses to create this. The main difference between my implementation is that I skip auxiliary variable s2.

Matlab Code for number of bodies set to 5:

% Adam Wermus

% July 1, 2016

% This program uses the Parker-Sochacki Method to solve the n-Body problem

% in 2 dimensions

% clear everything

close all;

clear all;

clc;

% Parameters

% num_of_bodies is the number of bodies

% G is the gravitational constant.

% Next are the masses for every “body”

num_of_bodies = 5;

G = 1;

m(1) = 1;

m(2) = 1;

m(3) = 1;

m(4) = 1;

m(5) = 1.5;

% Declare variables for Parker-Sochacki Implementation

order = 7; % order of Maclaurin series

t = 0.005; % time step

frames = 180; % number of iterations

power = 3; % power of a power series taking s to s^3

% Create arrays for all Maclaurin coefficients for position and velocity in

% the x and y directions

x_pos = zeros(num_of_bodies, order-1);

y_pos = zeros(num_of_bodies, order-1);

vx = zeros(num_of_bodies, order-1);

vy = zeros(num_of_bodies, order-1);

% Create arrays for the auxiliary variables s, s3, and a

s = zeros(num_of_bodies,num_of_bodies,order-1);

s3 = zeros(num_of_bodies,num_of_bodies,order-1);

a = zeros(num_of_bodies,num_of_bodies,order-1);

% initial mass 1 conditions

x_pos(1,1) = 0.55;

vx(1,1) = 0;

y_pos(1,1) = 0.27;

vy(1,1) = 0.23;

% initial mass 2 conditions

x_pos(2,1) = 0.30;

vx(2,1) = 0;

y_pos(2,1) = 0.27;

vy(2,1) = -2.42;

% initial mass 3 conditions

x_pos(3,1) = -0.30;

vx(3,1) = 0;

y_pos(3,1) = 0.27;

vy(3,1) = 2.42;

% initial mass 4 conditions

x_pos(4,1) = -0.55;

vx(4,1) = 0;

y_pos(4,1) = 0.27;

vy(4,1) = -0.23;

% initial mass 5 conditions

x_pos(5,1) = 0.0;

vx(5,1) = 0;

y_pos(5,1) = -0.72;

vy(5,1) = 0;

% initialize x position, y position, velocity x, velocity y for the bodies

x_pos_sum = zeros(num_of_bodies,frames);

y_pos_sum = zeros(num_of_bodies,frames);

vx_sum = zeros(num_of_bodies, frames);

vy_sum = zeros(num_of_bodies, frames);

% This loop gives the first frame values of position and velocity to the

% arrays

for n = 1:num_of_bodies

x_pos_sum(n,1) = x_pos(n,1);

y_pos_sum(n,1) = y_pos(n,1);

vx_sum(n,1) = vx(n,1);

vy_sum(n,1) = vy(n,1);

end

% This for loop updates the frames

for j = 2:frames

% Calculate Coefficients

for iterate = 1 : order-1

% initialize auxiliary variables

% This for loop guarantees that we won’t go to the same “body”

% This loop does half of the matrix calculation because it is

% symmetric(the distance between point 1 and point 2 is the same as

% the distance between point 2 and point 1

% We fill in the other parts in the next loop

for n = 1:num_of_bodies

for p = n+1:num_of_bodies

s(n,p,1) = 1./sqrt((x_pos(n,1) – x_pos(p,1)).^2 + (y_pos(n,1) – y_pos(p,1)).^2);

end

end

% This avoids doing multiple calculations over again. This is a

% symmetric matrix with diagonal zeros(actually the diagonals for

% this are infinity

for n = 1:num_of_bodies

for p = n+1:num_of_bodies

s(p,n,1) = s(n,p,1);

end

end

% Cauchy Product Arrays

cp_vx = zeros(num_of_bodies);

cp_vy = zeros(num_of_bodies);

cp_a = zeros(num_of_bodies, num_of_bodies);

cp_s = zeros(num_of_bodies, num_of_bodies); % cauchy product matrix for auxiliary variable s

cp_s3 = zeros(num_of_bodies, num_of_bodies); % cauchy product matrix for auxiliary variable s3

% Cauchy Products s3, v3, x3 first setup

if iterate == 1

% Computes initial condition of first half of matrix

for n = 1:num_of_bodies

for p = n+1:num_of_bodies

s3(n,p,1) = s(n,p,1) * s(n,p,1) * s(n,p,1);

end

end

% Fills in second half of matrix

for n = 1:num_of_bodies

for p = n+1:num_of_bodies

s3(p,n,1) = s3(n,p,1);

end

end

% For all other non-initial auxiliary variable coefficients

else

for k = 2:iterate

% Putting quick here avoids calculating this muliple times

quick = 1/(iterate – 1);

% Here we sum through the bodies to get our s3 term. This

% uses a power of a power series. We only do half of the

% matrix because it is symmetric. The other half is filled

% out later

for n = 1:num_of_bodies

for p = n+1:num_of_bodies

cp_s3(n,p) = cp_s3(n,p) + ((k-1) * power – (iterate – 1) + (k-1)) * (s(n, p, k) * s3(n, p, iterate + 1 – k));

end

end

end

% This part comes from the last chunk of the formula to take a

% power of a power series

for n = 1:num_of_bodies

for p = n+1:num_of_bodies

cp_s3(n,p) = cp_s3(n,p) * (quick / s(n,p,1)) ;

end

end

end

% The formula only applies to the non-initial condition

if iterate >= 2

for n = 1:num_of_bodies

for p = n+1:num_of_bodies

s3(n,p,iterate) = cp_s3(n,p);

end

end

% Fill in second half of s3 matrix

for n = 1:num_of_bodies

for p = n+1:num_of_bodies

s3(p,n,iterate) = s3(n,p,iterate);

end

end

end

% Cauchy Product for auxiliary variable a

% We cover only half of the matrix because it is symmetric

for k = 1:iterate

for n = 1:num_of_bodies

for p = n+1:num_of_bodies

cp_a(n,p) = cp_a(n,p) + ((x_pos(n,k) – x_pos(p,k)).*(vx(n,(iterate+1)-k)-vx(p,(iterate+1)-k)) + (y_pos(n,k) – y_pos(p,k)).*(vy(n,(iterate+1)-k)-vy(p,(iterate+1)-k)));

end

end

end

% Fills in second half of matrix

for n = 1:num_of_bodies

for p = n+1:num_of_bodies

cp_a(p,n) = cp_a(n,p);

end

end

% Fill in values into auxiliary variable a from cauchy product

% calculation above

for n = 1:num_of_bodies

for p = 1:num_of_bodies

a(n,p,iterate) = cp_a(n,p);

end

end

% Cauchy Product for auxiliary variable s

% We cover only half of the matrix because it is symmetric

for k = 1:iterate

for n = 1:num_of_bodies

for p = n+1:num_of_bodies

cp_s(n,p) = cp_s(n,p) + (s3(n,p,k)*a(n,p,(iterate+1)-k));

end

end

end

% Fills in second half of matrix s

for n = 1:num_of_bodies

for p = n+1:num_of_bodies

cp_s(p,n) = cp_s(n,p);

end

end

% Fill in values into auxiliary variable a from cauchy product

% calculation above

for n = 1:num_of_bodies

for p = 1:num_of_bodies

s(n,p,iterate + 1) = -cp_s(n,p)/iterate;

end

end

% Updates position for all bodies

for n = 1 : num_of_bodies

x_pos(n,iterate+1) = vx(n,iterate)/iterate;

y_pos(n,iterate+1) = vy(n,iterate)/iterate;

end

% Cauchy Products for Velocity

for k = 1:iterate

% We need to sum over all of the neighboring bodies. This is

% the Cauchy Product of the velocity

for n = 1:num_of_bodies

for p = 1:num_of_bodies

cp_vx(n) = cp_vx(n) + (G*m(p)*(x_pos(p,k)-x_pos(n,k))*s3(n,p,(iterate+1)-k));

cp_vy(n) = cp_vy(n) + (G*m(p)*(y_pos(p,k)-y_pos(n,k))*s3(n,p,(iterate+1)-k));

end

end

end

% We will now put the values from the caucy product into the

% velocity array

for n = 1 : num_of_bodies

vx(n,iterate+1) = cp_vx(n)/iterate;

vy(n,iterate+1) = cp_vy(n)/iterate;

end

end

% Setup sums for the bodies

x_temp_sum = zeros(num_of_bodies);

y_temp_sum = zeros(num_of_bodies);

vx_temp_sum = zeros(num_of_bodies);

vy_temp_sum = zeros(num_of_bodies);

% Compute Power series

for i = 1:order-1

for n = 1:num_of_bodies

x_temp_sum(n) = x_temp_sum(n) + (x_pos(n,i)*t^(i-1));

y_temp_sum(n) = y_temp_sum(n) + (y_pos(n,i)*t^(i-1));

vx_temp_sum(n) = vx_temp_sum(n) + (vx(n,i)*t^(i-1));

vy_temp_sum(n) = vy_temp_sum(n) + (vy(n,i)*t^(i-1));

end

end

% Put in power series to store in array

for n = 1:num_of_bodies

x_pos_sum(n,j) = x_temp_sum(n);

y_pos_sum(n,j) = y_temp_sum(n);

vx_sum(n,j) = vx_temp_sum(n);

vy_sum(n,j) = vy_temp_sum(n);

end

% Reset Initial Conditions

for n = 1:num_of_bodies

x_pos(n,1) = x_temp_sum(n);

y_pos(n,1) = y_temp_sum(n);

vx(n,1) = vx_temp_sum(n);

vy(n,1) = vy_temp_sum(n);

end

end

% Plot Results

% We know this is a 3-body version of n-body. This has to change as the

% number of bodies change.

% This creates 1-D Arrays for x and y of every “body”

x1 = x_pos_sum(1,:);

y1 = y_pos_sum(1,:);

x2 = x_pos_sum(2,:);

y2 = y_pos_sum(2,:);

x3 = x_pos_sum(3,:);

y3 = y_pos_sum(3,:);

x4 = x_pos_sum(4,:);

y4 = y_pos_sum(4,:);

x5 = x_pos_sum(5,:);

y5 = y_pos_sum(5,:);

% Plot results

plot(x1,y1,’c-‘, x2,y2, ‘k’, x3,y3,’r’, x4,y4, ‘b’, x5,y5,’g.’);

legend(‘mass 1’, ‘mass 2’, ‘mass 3’, ‘mass 4’, ‘mass 5’);

title(‘Parker-Sochacki N-Body Implementation where N = 5’);

Here is how this is broken down:

I have an x position matrix n by k where n is the particle and k is the order.

I have a y position matrix n by k where n is the particle and k is the order.

I have an x velocity matrix n by k where n is the particle and k is the order.

I have a y velocity matrix n by k where n is the particle and k is the order.

I have an auxiliary variable matrix for s,s3,and a that is n by n by k. This distinguishes the relationship between particles for order k.

I have matrices that are used to calculate the cauchy products for velocity, s, s3, and a. These matrices are n by n.

In the code, there are a lot of spots where I sum over all neighboring particles.

Here is an important note about the matrices for s, s3, and a. The relationship between particle 1 and particle 2 is the same as the relationship between particle 2 and particle 1. All 3 of these matrices are symmetric. I only need to compute half of those values(you can see this in the code where I fill half of one matrix and use that to fill in the other half). The diagonals for s, s3, and a were set to zero. The diagonals for s should be infinity because I would be dividing by zero.

Regarding s3, I still took the power of a power series as opposed to using cauchy products to compute s2, and then using cauchy products to compute s3. In a previous blog, I found that taking a power of a power series is faster. In all of the n-body Parker-Sochacki papers I have come across, I have not seen this type of implementation where s2 is skipped. I am not sure why.

When I implemented the 3-body replica in a previous blog, the simulation time took around 5 seconds as opposed to 1. The way I generalized this algorithm led to a greater computation cost.

Next I will show outputs for different initial conditions that I found online. I will show the image followed by the initial conditions in case you would be interested in trying this yourself.

This is the original 2-body problem that I showed a few months ago.

number of bodies = 2

G = 1

mass of particle 1 = 1

mass of particle 2 = 2

Order of Maclaurin series = 22

time step = 0.04167

number of frames = 10000

x_particle 1 = 1

velocity_x_particle 1 = 0.65

y_particle 1 = 0.5

velocity_y_particle 1 = 0.2

x_particle 2 = -1

velocity_x_particle 2 = -0.45

y_particle 2 = -0.3

velocity_y_particle 2 = 0.3

This is the original 3-body problem that I showed a few months ago.

number of bodies = 3

G = 1

mass of particle 1 = 0.6

mass of particle 2 = 0.1001

mass of particle 3 = 0.00101

Order of Maclaurin series = 30

time step = 0.001

number of frames = 2125

x_particle 1 = 0

velocity_x_particle 1 = 0.6

y_particle 1 = 1.95

velocity_y_particle 1 = 0.5

x_particle 2 = 0

velocity_x_particle 2 = 2.828

y_particle 2 = 2

velocity_y_particle 2 = 0

x_particle 3 = 0

velocity_x_particle 3 = 0

y_particle 3 = 2.1

velocity_y_particle 3 = 0

3-Body Problem where the red particle that cannot be seen is the center of mass. The blue and black particles are going around it

number of bodies = 3

G = 1

mass of particle 1 = 1

mass of particle 2 = 1

mass of particle 3 = 2

Order of Maclaurin series = 30

time step = 0.1

number of frames = 2000

x_particle 1 = -1

velocity_x_particle 1 = 0

y_particle 1 = 0

velocity_y_particle 1 = 1

x_particle 2 = 1

velocity_x_particle 2 = 0

y_particle 2 = 0

velocity_y_particle 2 = -1

x_particle 3 = 0

velocity_x_particle 3 = 0

y_particle 3 = 0

velocity_y_particle 3 = 0

This is the 3-body problem that simulates a figure 8. All 3 bodies are tracing the same path, which is why you only see one particle.

number of bodies = 3

G = 1

mass of particle 1 = 1

mass of particle 2 = 1

mass of particle 3 = 1

Order of Maclaurin series = 10

time step = 0.001

number of frames = 9000

x_particle 1 = 0.97000436

velocity_x_particle 1 = 0.466203685

y_particle 1 = -0.24308753

velocity_y_particle 1 = 0.43236573

x_particle 2 = -0.97000436

velocity_x_particle 2 = 0.466203685

y_particle 2 = 0.24308753

velocity_y_particle 2 = 0.43236573

x_particle 3 = 0

velocity_x_particle 3 = -0.93240737

y_particle 3 = 0

velocity_y_particle 3 = -0.86473146

This is a 3-body problem where I only modified the x position of particle 1 from the above simulation. I also increased the time step by a factor of 10. Everything else is the same resulting in a figure 8 that rotates.

number of bodies = 3

G = 1

mass of particle 1 = 1

mass of particle 2 = 1

mass of particle 3 = 1

Order of Maclaurin series = 10

time step = 0.01

number of frames = 9000

x_particle 1 = 0.95000436

velocity_x_particle 1 = 0.466203685

y_particle 1 = -0.24308753

velocity_y_particle 1 = 0.43236573

x_particle 2 = -0.97000436

velocity_x_particle 2 = 0.466203685

y_particle 2 = 0.24308753

velocity_y_particle 2 = 0.43236573

x_particle 3 = 0

velocity_x_particle 3 = -0.93240737

y_particle 3 = 0

velocity_y_particle 3 = -0.86473146

This is simply chaotic motion with 5 bodies

number of bodies = 5

G = 1

mass of particle 1 = 1

mass of particle 2 = 1

mass of particle 3 = 1

mass of particle 4 = 1

mass of particle 5 = 1

Order of Maclaurin series = 7

time step = 0.005

number of frames = 180

x_particle 1 = 0.55

velocity_x_particle 1 = 0

y_particle 1 = 0.27

velocity_y_particle 1 = 0.23

x_particle 2 = 0.3

velocity_x_particle 2 = 0

y_particle 2 = 0.27

velocity_y_particle 2 = -2.42

x_particle 3 = -0.30

velocity_x_particle 3 = 0

y_particle 3 = 0.27

velocity_y_particle 3 = 2.42

x_particle 4 = -0.55

velocity_x_particle 4 = 0

y_particle 4 = 0.27

velocity_y_particle 4 = -0.23

x_particle 5 = 0

velocity_x_particle 5 = 0

y_particle 5 = -0.72

velocity_y_particle 5 = 0

These are all of the examples I tried but you can play with the code to any number of particles you would like to get some interesting patterns!

I have a final note on the order and time step. If you look, some simulations had much smaller time steps than others. My smallest time step was 0.001 while my largest time step was 0.01. My lowest order was 7 while my highest order was 30. I was solving the same system of Ordinary Differential Equations. The only differences were in the initial conditions. If I were using other numerical methods like Euler, Verlet, Leap Frog, or Runge Kutta, I could only adjust the time step(I would have to rewrite Runge-Kutta to get methods like Runge-Kutta 2nd Order, Runge-Kutta 3rd Order, and Runge-Kutta 4th Order). If my order was too low with Parker-Sochacki, I only had to change that value and run the simulation again. With the Parker-Sochacki method, I was able to play with both time step and order just by changing what those values are set to.

I also came across the same “close encounters” benefit in the 5-body problem. If I lowered the order too much, the particles would all diverge. I could compensate by lowering the time step or raising the order of the Maclaurin series. I did this experimentally.

Take care and have a happy 4th of July!

Sincerely,

Adam Wermus