# Back to the 2-Body

In an earlier blog I applied the Parker-Sochacki method to the 2-Body gravity problem:

This was important because in games and film it is common to have a simulation where you have |r_1 – r_2| in the denominator. I’ve seen it in Smoothed Particle Hydrodynamics and Force Based Cloth. This blog showed that Parker-Sochacki could large time steps and have a high accuracy. Parker-Sochacki is a great fit when these Ordinary Differential Equations are stiff. For more details please refer to the blog a little ways down.

There is another application I want to show for this and it is in regards to particle simulations. Suppose one of those bodies is fixed so instead we have something like this:

x_0 and y_0 are fixed. I applied Parker-Sochacki to this simulation modifying initial conditions of x_0, y_0, gravity, mass, x(0), and y(0)(those are the positions of the masses that change where x_0 and y_0 are stationary not changing). This won’t sound scientific but I got cool different patterns that would be great particle effects in a film or video game. I have 10 examples followed by the code in case you would like to play with this in Matlab for yourself. I set the time step to 1/30 to simulate a video game atmosphere.

CODE:

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

% February 5, 2016

% This program uses the Parker-Sochacki Method to solve a 2D 2 body problem

% where 1 body is fixed

% clear everything

close all;

clear all;

clc;

% Parameters

G = 1;

m1 = 5;

% Stationary

x_0 = 0;

y_0 = 1;

% Declare variables for Parker-Sochacki Implementation

order = 10; % order for power series

%t = 0.04167./2; % time step

t = 0.03333333333;

frames = 2000; % number of iterations

% initial mass 1 conditions

%x_pos1(1) = 1;

%vx1(1) = 0.65;

%y_pos1(1) = 0.5;

%vy1(1) = 0.2;

x_pos1(1) = 1;

vx1(1) = -0.65;

y_pos1(1) = 1.5;

vy1(1) = 0.2;

% initial auxiliary variable conditions

r(1) = sqrt((x_0 – x_pos1(1)).^2 + (y_0 – y_pos1(1)).^2);

d_check(1) = sqrt((x_0 – x_pos1(1)).^2 + (y_0 – y_pos1(1)).^2);

s(1) = 1./r(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);

% Auxiliary Variable sums

r_sum = zeros(1, frames);

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

s_sum(1) = s(1);

r_sum(1) = r(1);

b(1) = (x_0 – x_pos1(1))*(vx1(1)) + (y_0 – y_pos1(1))*(vy1(1));

c(1) = x_0 * vx1(1) – x_pos1(1) * vx1(1) + y_0 * vy1(1) – y_pos1(1) * vy1(1);

% This for loop updates the frames

for j= 2:frames

% Calculate Coefficients

for iterate = 2 : order-1

% very first terms

s2(1) = s(1)*s(1);

s3(1) = s2(1)*s(1);

a(1) = (x_0 – x_pos1(1))*(vx1(1)) + (y_0 – y_pos1(1))*(vy1(1));

b(1) = (x_0 – x_pos1(1))*(vx1(1)) + (y_0 – y_pos1(1))*(vy1(1));

c(1) = x_0 * vx1(1) – x_pos1(1) * vx1(1) + y_0 * vy1(1) – y_pos1(1) * vy1(1);

% setup first terms

x_pos1(2) = vx1(1);

y_pos1(2) = vy1(1);

vx1(2) = G * m1 * x_0 – G * m1 * x_pos1(1) * s3(1);

vy1(2) = G * m1 * y_0 – G * m1 * y_pos1(1) * s3(1);

s(2) = -a(1) * s3(1);

r(2) = -a(1) * s(1);

% Setup cauchy product

% Velocity

cp_vx1 = 0;

cp_vy1 = 0;

% Auxiliary Cauchy Product

cp_s = 0;

cp_r = 0;

cp_s2 = 0;

cp_s3 = 0;

cp_a = 0;

cp_b = 0;

cp_b1 = 0;

cp_b2 = 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_0 – x_pos1(k)).*vx1(iterate+1-k) + (y_0 – y_pos1(k)).*vy1(iterate+1-k));

cp_b1 = cp_b1 + (x_pos1(k).*vx1(iterate+1-k));

cp_b2 = cp_b2 + (y_pos1(k).*vy1(iterate+1-k));

end

s3(iterate) = cp_s3; % s3

%a(iterate) = cp_a; % a

b(iterate) = x_0 * vx1(iterate) – cp_b1 + y_0 * vy1(iterate) – cp_b2;

a(iterate) = b(iterate); % a

% Position for mass 1

x_pos1(iterate+1) = vx1(iterate)/iterate;

y_pos1(iterate+1) = vy1(iterate)/iterate;

% Cauchy Products

% Compute Cauchy Products

for k =1:iterate

% Velocity for mass 1

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

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

% Auxiliary Cauchy Product

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

cp_r = cp_r + ((a(k)*s((iterate+1)-k))/iterate);

end

% % Put in cauchy product to coefficient

% Velocity for mass 1

vx1(iterate+1) = -cp_vx1;

vy1(iterate+1) = -cp_vy1;

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

r(iterate+1) = -cp_r;

end

%

% % Setup sums

x1_temp_sum = 0;

y1_temp_sum = 0;

vx1_temp_sum = 0;

vy1_temp_sum = 0;

s_temp_sum = 0;

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

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

r_temp_sum = r_temp_sum + (r(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;

s_sum(j) = s_temp_sum;

r_sum(j) = r_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;

s(1) = s_temp_sum;

r(1) = r_temp_sum;

d_check(j) = sqrt((x_0 – x_pos1(1)).^2 + (y_0 – y_pos1(1)).^2);

end

plot(x_pos1_sum,y_pos1_sum,x_0,y_0);

%plot(1:frames,s_sum, ‘b’, 1:frames,r_sum, ‘r’, 1:frames, d_check, ‘g’);

%legend(‘s variable’, ‘r variable’, ‘square root’);

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

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

Take care!