Greetings!

This post talks about applying the Parker-Sochacki method to the 2-body gravity problem. This means that there are two masses orbiting around each other. Implementation, code, and results are attached for a 2D and 3D version. I chose this simulation because it is the base case for the n-body gravity problem. I would also like to note that Dr. Rudmin, Dr. Pruett, Dr. Lacy, Dr. Ingham, and Dr. Herman have published papers on applying Parker-Sochacki to the n-body gravitational problem. They found a way to implement it and then make do it again in parallel. I applied Parker-Sochacki in the same way and will use the same substitutions that they did.

This was also a problem that was covered in the summer at BYU(the 2-body version).

I will start with the 2D version of the 2-body problem that was implemented into Matlab. Below are two images to show what it is supposed to look like. I chose a time step of 0.04167(close to 1/24) and, in Matlab, chose to set Parker-Sochacki to Order 20(note I said Matlab).

This first image goes for 1000 frames:

This second image goes for 10000 frames:

This shows you the pattern of what is supposed to happen. Two masses moved around each other as they traveled.

Now I would like to show what happened at lower orders, which is the same accuracy as the explicit methods like Forward Euler and Runge-Kutta. For this simulation, I am asking Matlab to update the position and velocity around 10,000 times at a time step that is typical for animation. Below are images over 10,000 frames as I increase the order starting at order 3 for Matlab(this means in C++ or Math it would be order 2).

The last few look similar and the accuracy is there. This simulation was not accurate until we achieved an order in the upper teens/low 20s. Only a higher order method(Parker-Sochacki) was able to handle such a large time step.

Okay another important issue I need to address. There are higher order methods out there. What good would it do if it is much slower? There will be a balance between what order is chosen and what time step is chosen. I can increase the order to 100 but it will be much slower. Here is a table of how long every simulation took in Matlab starting with order 3. I am documenting the “Self Time” check in Matlab. This is a mean value of 3 checks.

Order 3: 0.214 seconds

Order 4: 0.2813 seconds

Order 5: 0.3716 seconds

Order 6: 0.4753 seconds

Order 7: Â 0.588 seconds

Order 8: 0.7093 seconds

Order 9: 0.846 seconds

Order 10: 0.9926 seconds

Order 11: 1.1423 seconds

Order 12: 1.293 seconds

Order 13: 1.4873 seconds

Order 14: 1.6573 seconds

Order 15: 1.8473 seconds

Order 16: 2.0496 seconds

Order 17: 2.278 seconds

Order 18: 2.5083 seconds

Order 19: 2.7503 seconds

Order 20: 3.0026 seconds

Order 21: 3.27 seconds

Order 22: 3.5283 seconds

So yes, a higher order takes longer to compute. I don’t have access to Microsoft Excel here but I plugged this data into an online plot and this is the relationship of time vs order:

Out of curiosity, I checked how long it took for Order 100 and well:

Order 100: 56.8843 seconds

This is important because if I use a lower order to get the same accuracy as I did with the higher order Parker-Sochacki, I need to lower the time step and increase the number of frames to get the same simulation.

Remember all of this is in Matlab. The order numbers I talk about are for Matlab notation. They start at 1 not 0.

So when I lowered the order to 3, no time step could get an accurate result. At order 4, when I divided the time step(0.04167) by 7 and multiplied the number of frames by 7, I got the same accuracy as an order in the upper teens set by Parker-Sochacki. It took an average of 1.431 seconds. This took the amount of time between order 12 and order 13 and was more accurate.

But what if I use a higher order and decrease the time step? I set Parker-Sochacki to order 10, divided the time step by 2, and multiplied the number of frames by 2 and the simulation was completed in 1.0836 seconds(this is faster than order 4, which took 1.431 seconds).

That is a noticeable speed up by choosing to use Parker-Sochacki instead of the Runge-Kutta methods for this simulation. I could get a higher accuracy, take larger time steps, and improve the total simulation time by choosing the right order and time step.

I implemented a 3D version of this into Houdini where I multiplied the time step by 10 (so now it is 0.4167). Runge-Kutta 2nd Order crashed. Runge-Kutta 3rd Order was not accurate because the spheres moved around each other and then diverged away. Runge-Kutta 4th Order held the orbital motion. Parker-Sochacki was stable and accurate for the duration of the simulation.

2-Body version implemented with Runge-Kutta 2nd Order(it crashes):

2-Body version implemented with Parker-Sochacki at Order 10(it is stable):

Below is the Matlab Code for the 2D 2-body implementation:

————————————————————-

% Adam Wermus

% July 18, 2015

% This program uses the Parker-Sochacki Method to solve 6.1 2D

% clear everything

close all;

clear all;

clc;

% Parameters

G = 1;

m1 = 1;

m2 = 2;

% Declare variables for Parker-Sochacki Implementation

order = 22; % order for power series

t = 0.04167; % time step

frames = 10000; % number of iterations

% initial mass 1 conditions

x_pos1(1) = 1;

vx1(1) = 0.65;

y_pos1(1) = 0.5;

vy1(1) = 0.2;

% initial mass 2 conditions

x_pos2(1) = -1;

vx2(1) = -0.45;

y_pos2(1) = -0.3;

vy2(1) = 0.3;

% initial auxiliary variable conditions

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_pos2(1) – x_pos1(1))*(vx2(1) – vx1(1)) + (y_pos2(1) – y_pos1(1))*(vy2(1)-vy1(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);

% Auxiliary Variable sums

s_sum = zeros(1, frames);

s2_sum = zeros(1, frames);

s3_sum = zeros(1, frames);

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

s_sum(1) = s(1);

s2_sum(1) = s2(1);

s3_sum(1) = s3(1);

a_sum(1) = a(1);

% This for loop updates the frames

for j= 2:frames

% Setup cauchy product

cp_vx1 = 0;

cp_vy1 = 0;

cp_vx2 = 0;

cp_vy2 = 0;

cp_s = 0;

cp_s2 = 0;

cp_s3 = 0;

cp_a = 0;

% Calculate Coefficients

for iterate = 1 : order-1

% Setup cauchy product

% Velocity

cp_vx1 = 0;

cp_vy1 = 0;

cp_vx2 = 0;

cp_vy2 = 0;

% Auxiliary Cauchy Product

cp_s = 0;

cp_s2 = 0;

cp_s3 = 0;

cp_a = 0;

% Cauchy products

for k = 1:iterate

cp_s2 = cp_s2 + (s(k).*s(iterate+1-k));

end

% cauchy product iterates

s2(iterate) = cp_s2; % s2

for k = 1:iterate

cp_s3 = cp_s3 + (s2(k).*s(iterate+1-k));

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

end

s3(iterate) = cp_s3; % s3

a(iterate) = cp_a; % a

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

% Cauchy Products

% Compute Cauchy Products

for k =1:iterate

% Velocity for mass 1

cp_vx1 = cp_vx1 + ((-G*m2*(x_pos1(k)-x_pos2(k))*s3((iterate+1)-k))/iterate);

cp_vy1 = cp_vy1 + ((-G*m2*(y_pos1(k)-y_pos2(k))*s3((iterate+1)-k))/iterate);

% Velocity for mass 2

cp_vx2 = cp_vx2 + ((-G*m1*(x_pos2(k)-x_pos1(k))*s3((iterate+1)-k))/iterate);

cp_vy2 = cp_vy2 + ((-G*m1*(y_pos2(k)-y_pos1(k))*s3((iterate+1)-k))/iterate);

% Auxiliary Cauchy Product

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

cp_s2 = cp_s2 + ((s(k)*s((iterate+1)-k)));

cp_s3 = cp_s3 + ((s2(k)*s((iterate+1)-k)));

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

end

% % Put in cauchy product to coefficient

% Velocity for mass 1

vx1(iterate+1) = cp_vx1;

vy1(iterate+1) = cp_vy1;

% Velocity for mass 1

vx2(iterate+1) = cp_vx2;

vy2(iterate+1) = cp_vy2;

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

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;

s_temp_sum = 0;

s2_temp_sum = 0;

s3_temp_sum = 0;

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

s_temp_sum = s_temp_sum + (s(i)*t^(i-1));

s2_temp_sum = s2_temp_sum + (s2(i)*t^(i-1));

s3_temp_sum = s3_temp_sum + (s3(i)*t^(i-1));

a_temp_sum = a_temp_sum + (a(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;

s_sum(j) = s_temp_sum;

s2_sum(j) = s2_temp_sum;

s3_sum(j) = s3_temp_sum;

a_sum(j) = a_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;

s(1) = s_temp_sum;

s2(1) = s2_temp_sum;

s3(1) = s3_temp_sum;

a(1) = a_temp_sum;

end

plot(x_pos1_sum,y_pos1_sum,x_pos2_sum,y_pos2_sum);

title(‘Order 22 at a Time Step of 0.04167’);

————————————————————-

My last note is to address why this entire example would matter in Computer Animation. First off, this can be expanded into an n-body orbital problem, which is something I am implementing into Houdini at the moment. It could create some interesting effects. Another reason is because of what the coupled Ordinary Differential Equations are:

The denominator where yo have the difference in position, |r1-r2|, is the auxiliary variable to use Parker-Sochacki. That setup is common in physical simulations. I have seen it in Force based Cloth and Smoothed Particle Hydrodynamics. The structure is similar.

I applied an auxiliary variable to |r1-r2| and created a system of Ordinary Differential Equations that could produce the following recurrence relationship:

On the next episode of this blog, I will address a Parker-Sochacki recurrence relationship setup issue that has appeared twice now. The first time it came up with the Projectile with Gravity and Drag and it reappears here. I am talking about when the recurrence relationship is to compute i instead of i+1.

Thank you for going through this! So far this should be the most promising indication to implement Parker-Sochacki into physical simulations over other explicit methods.