This blog continues from the previous blog where I applied the Parker-Sochacki method to the 3-body gravity problem preparing for the n-body problem.

I answered the questions I posed in the last blog: Is it better to update the auxiliary variables as a Maclaurin series or do the initial condition calculation again? What is faster when taking a Maclaurin series to the 3rd power-Doing cauchy products twice or the power of a Maclaurin series formula? I also cut auxiliary variables, which saved time.

First I would like to point out that I am not changing the order or the time step. This is just optimizing the previous code. This is set to 11th order Parker-Sochacki(this means Matlab will show 12) at a time step of 0.001.

I timed the simulation 3 times and then took the mean value. I’ll show everything I did that lowered the time it took to run the simulation. I want to point out that this is for 3 bodies only. The improvements in time will be more significant as more bodies are part of the simulation.

On the Windows computer I am using at the University of Utah, the mean time from the previous code at 11th order is 1.1843 seconds. This is larger if I increase the order.

The first test was to see what is faster:using a Maclaurin series evaluation for the auxiliary variables s,u,v,w,x,y or using the initial conditions again. At order 11, using the square root auxiliary variables lowered the mean time to 1.1683 seconds. The accuracy was the same. In deciding what is more optimal, this depends on the order. If a higher order method is used, then most likely it will be better to reset the variables by calculating them the same way as when they were being set as initial conditions.

The next step was cutting auxiliary variables. I mentioned in the last blog that s and u are equal, v and w are equal, and x and y are equal. I cut u,w,y, which also cuts u2,w2,y2,u3,w3, and y3. This lowered the mean time to 1.0633 seconds.

The last thing I checked was taking the power of a power series that I showed towards the end of the last blog. This had to be done differently in Matlab because arrays start at 1. Taking the power of a power series took me straight to s3,v3,and x3 without needing s2,v2, and x2. This cut cauchy products. This lowered the mean time to 1.0516 seconds. From this, I would say generally taking the power of a power series (power 3) is faster than doing cauchy products to get s2,v2,x2 and then again to get s3,v3,x3.

Below is the latest code implemented in Matlab:

————————————————————————–

% Adam Wermus

% 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;

% extra for power of a power series

power = 3;

% Declare variables for Parker-Sochacki Implementation

order = 12; % 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));

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));

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));

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);

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);

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);

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);

v_sum(1) = v(1);

x_sum(1) = x(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;

cp_b = 0;

% Auxiliary Variables Cauchy Product 1-3

cp_v = 0;

cp_v2 = 0;

cp_v3 = 0;

cp_c = 0;

cp_d = 0;

% Auxiliary Variables Cauchy Product 2-3

cp_x = 0;

cp_x2 = 0;

cp_x3 = 0;

cp_e = 0;

cp_f = 0;

% Cauchy Products s3, v3, x3 first setup

if iterate == 1

s3(iterate) = s(iterate) * s(iterate) * s(iterate);

v3(iterate) = v(iterate) * v(iterate) * v(iterate);

x3(iterate) = x(iterate) * x(iterate) * x(iterate);

else

for k = 2:iterate

% First Part of Power of a Power Series Implementation in Matlab

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

cp_v3 = cp_v3 + ((k-1) * power – (iterate – 1) + (k-1)) * (v(k)*v3(iterate + 1 – k));

cp_x3 = cp_x3 + ((k-1) * power – (iterate – 1) + (k-1)) * (x(k)*x3(iterate + 1 – k));

end

% Second Part of Power of a Power Series Implementation in Matlab

quick = 1/(iterate – 1);

cp_s3 = cp_s3 * (quick / s(1));

cp_v3 = cp_v3 * (quick / v(1));

cp_x3 = cp_x3 * (quick / x(1));

end

if iterate >= 2

% Cauchy Product update s3, v3, x3

s3(iterate) = cp_s3; % s3

v3(iterate) = cp_v3; % v3

x3(iterate) = cp_x3; % x3

end

% 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, v, x

for k =1:iterate

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

cp_v = cp_v + (v3(k)*c((iterate+1)-k));

cp_x = cp_x + (x3(k)*e((iterate+1)-k));

end

% Auxiliary Variable separation update

s(iterate+1) = -cp_s/iterate; % s

v(iterate+1) = -cp_v/iterate; % v

x(iterate+1) = -cp_x/iterate; % x

% 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))*s3((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))*s3((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))*v3((iterate+1)-k));

cp_vx3 = cp_vx3 – (G*m2*(x_pos3(k)-x_pos2(k))*x3((iterate+1)-k));

% Velocity y for mass 3

cp_vy3 = cp_vy3 – (G*m1*(y_pos3(k)-y_pos1(k))*v3((iterate+1)-k));

cp_vy3 = cp_vy3 – (G*m2*(y_pos3(k)-y_pos2(k))*x3((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;

% 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));

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;

%

%

% % 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) = 1./sqrt((x_pos1(1) – x_pos2(1)).^2 + (y_pos1(1) – y_pos2(1)).^2);

v(1) = 1./sqrt((x_pos1(1) – x_pos3(1)).^2 + (y_pos1(1) – y_pos3(1)).^2);

x(1) = 1./sqrt((x_pos2(1) – x_pos3(1)).^2 + (y_pos2(1) – y_pos3(1)).^2);

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’);

Here is the attached image(this is the same as the previous blog):

This is a slightly more polished implementation than the previous version. Take care!