Parker-Sochacki 2-Body Gravity Problem


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:

PS_2_body 1000

This second image goes for 10000 frames:

PS_2_body 10000

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

Order 3 2-body

Order 4 2-body

Order 5 2-body

Order 6 2-body

Order 7 2-body.PNG

Order 8 2-body.PNG

Order 9 2-body.PNG

Order 10 2-body.PNG

Order 11 2-body.PNG

Order 12 2-body.PNG

Order 13 2-body.PNG

Order 14 2-body.PNG

Order 15 2-body.PNG

Order 16 2-body.PNG

Order 17 2-body.PNG

Order 18 2-body.PNG

Order 19 2-body.PNG

Order 20 2-body.PNG

Order 21 2-body.PNG

Order 22 2-body.PNG

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:

plot of time vs order.PNG

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;

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



Leave a Reply

Fill in your details below or click an icon to log in: Logo

You are commenting using your account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s