FIRST NOTE: There has been a very helpful website called myphysicslab(URL below) that has a lot of simulations I have used as a reference so that I can solve the ODEs using the Parker-Sochacki method:

Now onto the blog:

Greetings! This is the 4th part of the springs simulations. Before Springs Part 1, I wrote a blog in January about comparing Parker-Sochacki at order 4 and Runge-Kutta 4th order for a single spring in one dimension. I compared both numerical methods in terms of speed and accuracy in C++.

In Springs parts 1, 2, and 3 I showed how to apply the Parker-Sochacki method to the Ordinary Differential Equations for that specific spring simulation, showed code written in Matlab, and an animation implemented in Houdini. Part 1 applied to a single spring moving in 1 dimension. Part 2 applied to a single spring moving in 2 dimensions. Part 3 applied to a double spring moving in 1 dimension. In part 3, I also showed that a lower order Maclaurin series caused the “bobs” to “blow up” moving farther away from each other at the larger time step. Now part 4 will apply to a double spring moving in 2 dimensions.

I’ll start with a rendered video to show the simulation I am referring to. This was implemented and rendered in Houdini. I would like to note that, in this video, damping was set to 0 for both springs so they will not go to rest. Damping is turned on for the later videos.

For this blog, I will show the math on how I applied the Parker-Sochacki method. Then I will show code implemented in Matlab, and more videos rendered in Houdini.

We start with 8 variables, x_1, y_1, u_1, and v_1, x_2, y_2, u_2, and v_2. x_1 and y_1 refer to the x and y position of bob 1. u_1 and v_1 refer to the x velocity and y velocity respectively for bob 1. x_2 and y_2 refer to the x and y position of bob 2. u_2 and v_2 refer to the x velocity and y velocity respectively for bob 2. Our goal is to represent all 8 variables as a polynomial up to any order desired that is the Maclaurin series:

To do this, we will create extra variables called auxiliary variables. They will be used to develop a recursion relationship making it possible to compute as many coefficient terms up to order n as we want.

I will use some of the same notation that they used for the double spring.

First I will define two equations, L_1 and L_2:

where x_0 and y_0 are constants that represent the anchor point for where the top spring is attached.

I will now show the 8 First Order ODEs that describe the motion of the double spring in 2D and then go into how using the auxiliary variables work.

where k_1, k_2, m_1, m_2, r_1, r_2, b_1, and b_2, and g are constants. k_1 and k_2 are spring constants for spring 1 and spring 2 respectively. m_1 and m_2 are the masses of bob 1 and bob 2 respectively. r_1 and r_2 are the rest lengths for spring 1 and spring 2 respectively. b_1 and b_2 are the damping constants for spring 1 and spring 2 respectively. g is the gravity constant set to -9.81 m/s^2.

We introduce two auxiliary variables, a, and b. Let a = 1 / L_1 and b = 1/ L_2. They will also be converted into two polynomials that are also the Maclaurin series. After some algebra we can modify the velocity ODEs. We also need to take the derivatives of a and b. Before I do that, I need to introduce two new variables , c, and d, that will be part of a’ and b’ but I will not need to take their derivatives.

Now I will take the derivatives of a and b:

With this, we we now have a system of First Order Ordinary Differential Equations that can be used to create a system of Maclaurin series for the 8 variables, a, and b. For convenience, I will combine constants.

Let:

alpha = k_1 / m_1

beta = b_1 / m_1

gamma = k_2 / m_1

e = k_2 / m_2

f = b_2 / m_2

I will now show all 8 ODEs together with the algebra modifications to the velocity ODEs.

This is what we will use to create a Maclaurin coefficients up to order n. We start with initial positions and velocities. The initial values for a and b come from their equations where a = 1 / L_1 and b = 1 / L _2. Then using formulas from Power series, we create can always calculate coefficient i + 1 if we know coefficient i, which we always do.

2 things: the first coefficient for velocities are calculated differently than i + 1 when i is greater than or equal to 1 because of constants. I expand the terms from what is shown above. Another note is that we are going to use cauchy products where I will not multiply more than 2 variables at a time. To do this, I introduce more variables: b2, b3, a2, and a3. This is why c and d were created. You will see that I compute the ith term for b2, b3, a2, a3, c, and d. I will display computing the i+1 terms in order of position, velocity, auxiliary variables a and b, then a2, a3, b2, b3, c, and d:

Position:

Velocity:

Variables a and b:

The ith terms a2, a3, b2, b3, c, and d:

We are able to compute the 0th term for a2, a3, b2, b3, c, and d not needing their derivatives. We compute these coefficients to create a Maclaurin series for the 10 variables, 4 positions, 4 velocities, auxiliary variable a, and auxiliary variable b:

Once this Maclaurin series is calculated, they become the new initial conditions for the positions, velocities, and auxiliary variables. This continues until the simulation is done.

I want to make a special note about a and b. There are two choices on updating them. One option is presented here where the equations a = 1 / L_1 and b = 1 / L_2 are used once and for the rest of the simulation they are updated with a Maclaurin series. The other option is to use a = 1 / L_1 and b = 1 / L_2 for every frame. In a previous blog, I presented optimization options to the 3-body problem in preparing to become the n-body problem. The faster option depends on the order of the Maclaurin series.

I will now show the code implemented in Matlab. Unlike the previous springs, there is no “pretty picture” pattern.

% Adam Wermus

% July 13, 2016

% This solves the equations of motion for the 2-D Double Spring using the

% Parker-Sochacki method

% clear everything

clear all;

clc;

% Parameters

r_1 = 1.0; % rest length

r_2 = 1.0; % rest length

x_0 = 0; % anchor point x

y_0 = 2; % anchor point y

g = -9.81; % gravity

k_1 = 6.0; % spring stiffness

k_2 = 6.0; % spring stiffness

m_1 = 0.5; % mass of bob 1

m_2 = 0.5; % mass of bob 2

damping_1 = 0; % damping constant

damping_2 = 0; % damping constant

alpha = k_1 / m_1;

beta = damping_1 / m_1;

gamma = k_2 / m_1;

e = k_2 / m_2;

f = damping_2 / m_2;

% Parker Sochacki setup

order = 10;% order of Maclaurin series

h = 0.04167; % time step

frames = 5000; % number of frames;

% Setup Arrays for Maclaurin coefficients

x_1 = zeros(1, order);

y_1 = zeros(1, order);

x_2 = zeros(1, order);

y_2 = zeros(1, order);

u_1 = zeros(1, order);

v_1 = zeros(1, order);

u_2 = zeros(1, order);

v_2 = zeros(1, order);

a = zeros(1, order);

a2 = zeros(1, order);

a3 = zeros(1, order);

c = zeros(1, order);

b = zeros(1, order);

b2 = zeros(1, order);

b3 = zeros(1, order);

d = zeros(1, order);

% Setup Arrays that store Position and Velocity values

x_1_sum = zeros(1, frames);

y_1_sum = zeros(1, frames);

x_2_sum = zeros(1, frames);

y_2_sum = zeros(1, frames);

u_1_sum = zeros(1, order);

v_1_sum = zeros(1, order);

u_2_sum = zeros(1, order);

v_2_sum = zeros(1, order);

a_sum = zeros(1, frames);

b_sum = zeros(1, frames);

% Initial Conditions

x_1(1) = 0.5;

y_1(1) = -1.13333;

x_2(1) = 0;

y_2(1) = -2.45;

u_1(1) = 0;

v_1(1) = 0;

u_2(1) = 0;

v_2(1) = 0;

a(1) = 1./((x_1(1) – x_0).^2 + (y_1(1) – y_0).^2).^(1./2);

b(1) = 1./((x_2(1) – x_1(1)).^2 + (y_2(1) – y_1(1)).^2).^(1./2);

% Initial arrays

x_1_sum(1) = x_1(1);

y_1_sum(1) = y_1(1);

x_2_sum(1) = x_2(1);

y_2_sum(1) = y_2(1);

u_1_sum(1) = u_1(1);

v_1_sum(1) = v_1(1);

u_2_sum(1) = u_2(1);

v_2_sum(1) = v_2(1);

a_sum(1) = a(1);

b_sum(1) = b(1);

% Iterate through every frame

for j = 2:frames

% 1st

a2(1) = a(1).*a(1);

a3(1) = a2(1).*a(1);

b2(1) = b(1).*b(1);

b3(1) = b2(1).*b(1);

c(1) = ((x_1(1) – x_0).*u_1(1)) + (y_1(1) – y_0).*v_1(1);

d(1) = (x_2(1) – x_1(1)).*(u_2(1) – u_1(1)) + (y_2(1) – y_1(1)).*(v_2(1) – v_1(1));

% Second setup

x_1(2) = u_1(1);

y_1(2) = v_1(1);

x_2(2) = u_2(1);

y_2(2) = v_2(1);

% SPLIT UP u1 and v1 calculation

u_1(2) = -alpha.*x_1(1) + alpha.*x_0 + alpha.*r_1.*x_1(1).*a(1) – alpha.*r_1.*x_0.*a(1) – beta.*u_1(1) + gamma.*x_2(1) – gamma.*x_1(1) – gamma.*r_2.*x_2(1).*b(1) + gamma.*r_2.*x_1(1).*b(1);

v_1(2) = -alpha.*y_1(1);

v_1(2) = v_1(2) + alpha.*y_0;

v_1(2) = v_1(2) + alpha.*r_1.*y_1(1).*a(1);

v_1(2) = v_1(2) – alpha.*r_1.*y_0.*a(1);

v_1(2) = v_1(2) – beta.*v_1(1);

v_1(2) = v_1(2) + gamma.*y_2(1);

v_1(2) = v_1(2) – gamma.*y_1(1);

v_1(2) = v_1(2) – gamma.*r_2.*y_2(1).*b(1);

v_1(2) = v_1(2) + gamma.*r_2.*y_1(1).*b(1);

v_1(2) = v_1(2) + g;

u_2(2) = -e.*(x_2(1) – x_1(1)) + e.*r_2.*(x_2(1) – x_1(1)).*b(1) – f.*u_2(1);

v_2(2) = -e.*(y_2(1) – y_1(1)) + e.*r_2.*(y_2(1) – y_1(1)).*b(1) – f.*v_2(1) + g;

a(2) = -a3(1).*c(1);

b(2) = -b3(1).*d(1);

% Maclaurin coefficients

for i = 2:order-1

% cauchy product for a2 and b2

cp_a2 = 0;

cp_b2 = 0;

for k = 1:i

cp_a2 = cp_a2 + (a(k).*a(i+1-k));

cp_b2 = cp_b2 + (b(k).*b(i+1-k));

end

a2(i) = cp_a2;

b2(i) = cp_b2;

% cauchy product for b3

cp_a3 = 0;

cp_b3 = 0;

for k = 1:i

cp_a3 = cp_a3 + (a2(k).*a(i+1-k));

cp_b3 = cp_b3 + (b2(k).*b(i+1-k));

end

a3(i) = cp_a3;

b3(i) = cp_b3;

% cauchy product for c

cp_c = 0;

cp_d = 0;

for k = 1:i

cp_c = cp_c + ((x_1(k)-x_0).*u_1(i+1-k));

cp_c = cp_c + ((y_1(k)-y_0).*v_1(i+1-k));

cp_d = cp_d + ((x_2(k)-x_1(k)).*(u_2(i+1-k)-u_1(i+1-k)));

cp_d = cp_d + ((y_2(k)-y_1(k)).*(v_2(i+1-k)-v_1(i+1-k)));

end

c(i) = cp_c;

d(i) = cp_d;

% Update x and y position

x_1(i+1) = u_1(i) / i;

y_1(i+1) = v_1(i) / i;

x_2(i+1) = u_2(i) / i;

y_2(i+1) = v_2(i) / i;

% Velocity update u1 and v1

cp_u_1_1 = 0;

cp_v_1_1 = 0;

% First Part u1 and v1

for k = 1:i

cp_u_1_1 = cp_u_1_1 + (x_1(k).*a(i+1-k));

cp_v_1_1 = cp_v_1_1 + (y_1(k).*a(i+1-k));

end

cp_u_1_1 = alpha.*r_1.*cp_u_1_1;

cp_v_1_1 = alpha.*r_1.*cp_v_1_1;

cp_u_1_2 = 0;

cp_v_1_2 = 0;

% Second Part u1 and v1

for k = 1:i

cp_u_1_2 = cp_u_1_2 + (x_2(k).*b(i+1-k));

cp_v_1_2 = cp_v_1_2 + (y_2(k).*b(i+1-k));

end

cp_u_1_2 = -gamma.*r_2.*cp_u_1_2;

cp_v_1_2 = -gamma.*r_2.*cp_v_1_2;

cp_u_1_3 = 0;

cp_v_1_3 = 0;

% Third Part u1 and v1

for k = 1:i

cp_u_1_3 = cp_u_1_3 + (x_1(k).*b(i+1-k));

cp_v_1_3 = cp_v_1_3 + (y_1(k).*b(i+1-k));

end

cp_u_1_3 = gamma.*r_2.*cp_u_1_3;

cp_v_1_3 = gamma.*r_2.*cp_v_1_3;

%final push of u1 and v1

cp_u_1 = 0;

cp_v_1 = 0;

% NON CAUCHY PRODUCT PUSH

cp_u_1 = cp_u_1 + cp_u_1_1 + cp_u_1_2 + cp_u_1_3;

cp_v_1 = cp_v_1 + cp_v_1_1 + cp_v_1_2 + cp_v_1_3;

cp_u_1 = cp_u_1 – alpha.*x_1(i);

cp_v_1 = cp_v_1 – alpha.*y_1(i);

cp_u_1 = cp_u_1 – alpha.*r_1.*x_0.*a(i);

cp_v_1 = cp_v_1 – alpha.*r_1.*y_0.*a(i);

cp_u_1 = cp_u_1 – beta.*u_1(i);

cp_v_1 = cp_v_1 – beta.*v_1(i);

cp_u_1 = cp_u_1 + gamma.*(x_2(i) – x_1(i));

cp_v_1 = cp_v_1 + gamma.*(y_2(i) – y_1(i));

% Velocity update u1 and v1

u_1(i+1) = cp_u_1 / i;

v_1(i+1) = cp_v_1 / i;

% Velocity update u2 and v2

cp_u_2 = 0;

cp_v_2 = 0;

% First Part u2 and v2

for k = 1:i

cp_u_2 = cp_u_2 + (x_2(k)-x_1(k)).*b(i+1-k);

cp_v_2 = cp_v_2 + (y_2(k)-y_1(k)).*b(i+1-k);

end

cp_u_2 = e.*r_2.*cp_u_2;

cp_v_2 = e.*r_2.*cp_v_2;

% NON CAUCHY PRODUCT FOR u2

cp_u_2 = cp_u_2 – e.*(x_2(k)-x_1(k));

cp_v_2 = cp_v_2 – e.*(y_2(k)-y_1(k));

cp_u_2 = cp_u_2 – f.*u_2(i);

cp_v_2 = cp_v_2 – f.*v_2(i);

% Velocity Update u2 and v2

u_2(i+1) = cp_u_2 / i;

v_2(i+1) = cp_v_2 / i;

% cauchy product for a and b

cp_a = 0;

cp_b = 0;

for k = 1:i

cp_a = cp_a + (a3(k).*c(i+1-k));

cp_b = cp_b + (b3(k).*d(i+1-k));

end

a(i+1) = -cp_a / i;

b(i+1) = -cp_b / i;

end

% Temporary sums

x_1_sum_temp = 0;

y_1_sum_temp = 0;

x_2_sum_temp = 0;

y_2_sum_temp = 0;

u_1_sum_temp = 0;

v_1_sum_temp = 0;

u_2_sum_temp = 0;

v_2_sum_temp = 0;

a_sum_temp = 0;

b_sum_temp = 0;

% Compute Maclaurin series

for i = 1:order-1

x_1_sum_temp = x_1_sum_temp + (x_1(i)*h^(i-1));

y_1_sum_temp = y_1_sum_temp + (y_1(i)*h^(i-1));

x_2_sum_temp = x_2_sum_temp + (x_2(i)*h^(i-1));

y_2_sum_temp = y_2_sum_temp + (y_2(i)*h^(i-1));

u_1_sum_temp = u_1_sum_temp + (u_1(i)*h^(i-1));

v_1_sum_temp = v_1_sum_temp + (v_1(i)*h^(i-1));

u_2_sum_temp = u_2_sum_temp + (u_2(i)*h^(i-1));

v_2_sum_temp = v_2_sum_temp + (v_2(i)*h^(i-1));

a_sum_temp = a_sum_temp + (a(i)*h^(i-1));

b_sum_temp = b_sum_temp + (b(i)*h^(i-1));

end

% Store Maclaurin series into sum array

x_1_sum(j) = x_1_sum_temp;

y_1_sum(j) = y_1_sum_temp;

x_2_sum(j) = x_2_sum_temp;

y_2_sum(j) = y_2_sum_temp;

u_1_sum(j) = u_1_sum_temp;

v_1_sum(j) = v_1_sum_temp;

u_2_sum(j) = u_2_sum_temp;

v_2_sum(j) = v_2_sum_temp;

a_sum(j) = a_sum_temp;

b_sum(j) = b_sum_temp;

% Reset initial conditions

x_1(1) = x_1_sum(j);

y_1(1) = y_1_sum(j);

x_2(1) = x_2_sum(j);

y_2(1) = y_2_sum(j);

u_1(1) = u_1_sum(j);

v_1(1) = v_1_sum(j);

u_2(1) = u_2_sum(j);

v_2(1) = v_2_sum(j);

a(1) = a_sum(j);

b(1) = b_sum(j);

end

% Plot Results

%plot(x_1_sum,y_1_sum, ‘k-‘);

%title(‘Parker-Sochacki First BOB Double Spring 2D’)

plot(x_2_sum,y_2_sum, ‘k-‘);

title(‘Parker-Sochacki Second BOB Double Spring 2D’)

xlabel(‘x position’);

ylabel(‘y position’);

legend(‘Motion of a 2D Bob’);

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

The first plot is the motion for the first bob and the second plot is the motion for the second bob. Damping was set to 0. The order was set to 9(In Matlab this is 10 because arrays start at 1). The time step is 1/24. This went on for 5000 frames.

This simulation was also implemented in Houdini using the solver node on the bobs and a combination of the static solver and wire solver for the springs. In the first video, the damping was set to 0 so that the bobs would continue moving.

In this video, damping is set to 1 for both springs so they quickly move to rest. Also, the bottom bob’s velocity in x and y are set to 0. A greater value is set to the top bob in the x direction so that it starts moving to the right pulling the second bob with it. Quickly, damping pulls the bobs to the center where there is still a little up and down oscillation at the end. The time step is set to 0.054171.

In the other spring simulations, I set the time step to 1/12. While this simulation still moved accurately, it looked to fast. For these videos, I was shooting for a more natural look.

In this video, damping is set to 0.5 for the top spring and 0 to the bottom spring. The top bob is set at the origin with 0 velocity in both directions. The bottom bob is given a velocity of 8 in the x direction so it will take off to the right before it is pulled back. The velocity is 0 in the y direction:

In this video, damping is applied to both springs so that that the simulation comes to rest:

For all of these simulations, the mass for bob 1 and bob 2 were fixed at 0.5.

Thank you for taking the time to go through this! This has shown that Parker-Sochacki can be applied to 4 different types of spring simulations at higher time steps. This hasn’t covered collisions yet.