Greetings! In Springs Part 1, I applied the Parker-Sochacki method to solving the equations of motion for a single 1D spring. I also referenced a blog I wrote in January where I compared Parker-Sochacki 4th Order and Runge-Kutta 4th Order to the single spring in terms of speed and accuracy in C++. In Springs Part 2, I applied the Parker-Sochacki method to solving the equations of motion for a single 2D spring. For both parts, I showed the Ordinary Differential Equations, I showed the process of using the Parker-Sochacki method to convert a system of First Order Ordinary Differential Equations into a Maclaurin series up to any order desired, I showed code in Matlab, and I showed rendered videos of implementing these simulations in Houdini.

For Springs Part 3, we are applying the Parker-Sochacki method to a double spring in one dimension.

We’ll start with the rendered video in Houdini. This involves 2 “bobs” with 3 springs. Two springs are attached to the wall and “bob” and one other spring is attached to both bobs.

Here are the First Order Ordinary Differential Equations that make up the motion for the double spring in one dimension:

k_1, k_2, m_1, m_2, r_1, r_2, w_1, and w_2 are constants. k_1 and k_2 are spring constants. m_1 is the mass of bob 1 and m_2 is the mass of bob 2. r_1 is the rest length for bob 1. r_2 is the rest length for bob 2. w_1 is the width of bob 1. w_2 is the width of bob 2. x_1 is the position of bob 1. x_2 is the position of bob 2. u_1 is the velocity of bob 1. u_2 is the velocity of bob 2. This is all combined into the 4 Ordinary Differential Equations but if we wanted to just know the stretch of the spring, we would take x_1 – r_1 for the spring connecting the wall and bob 1. We would take x_2 – x_1 – w_1 – r_2 for the spring connecting bob 1 and bob 2.

Unlike Springs part 2 where we solved for the single spring in 2D, this simulation does not require creating extra equations nor does it require the use of cauchy products.

We start with an initial set of conditions for x_1, x_2, u_1, and u_2. I will also combine some constants together. Let

alpha = k_1 / m_1;

beta = k_2 / mass_1;

gamma = k_2 / mass_2;

The first subscript after x and u is the number of the bob. The only options are 1 and 2. This means that for x_2, it is the i + 1 term for bob 2, not multiplying i by 2 and adding 1. Also, like the previous blog, The first velocity terms are computed differently because there are constants being added.

Once all of the Maclaurin coefficients are computed, we have a polynomial for x_1, x_2, u_1, and u_2 that are the Maclaurin series or Taylor series centered at 0.

Once all 4 Maclaurin series are computed, they become the initial conditions for the next frame that is calculated. I will now show the code implemented in Matlab. I chose to set the order to 9(in matlab it will say 10 because arrays start at 1). The time step is set to 1/24.

% Adam Wermus

% July 9, 2016

% This solves the equations of motion for the Double Spring in 1D

% Parker-Sochacki method

% clear everything

clear all;

clc;

% constants

mass_1 = 1;

mass_2 = 1;

k_1 = 6; % spring stiffness 1

k_2 = 6; % spring stiffness 2

r_1 = 3; % length of spring

r_2 = 3; % length of spring

w_1 = 1; % width of block

w_2 = 1; % width of block

% simplification

alpha = k_1 / mass_1;

beta = k_2 / mass_1;

gamma = k_2 / mass_2;

% Parker-Sochacki Setup

order = 10; % order of Maclaurin series

frames = 2000; % number of frames

h = 0.04167; % time step

% Setup Arrays for Maclaurin coefficients

x_1 = zeros(1, order);

x_2 = zeros(1, order);

v_1 = zeros(1, order);

v_2 = zeros(1, order);

% Setup Arrays that store Position and Velocity values

x_1_sum = zeros(1, frames);

x_2_sum = zeros(1, frames);

v_1_sum = zeros(1, frames);

v_2_sum = zeros(1, frames);

% Initial Conditions

x_1(1) = 3.2;

x_2(1) = 6.4;

v_1(1) = 0.0;

v_2(1) = 0.0;

% Initial setup for sums

x_1_sum(1) = x_1(1);

x_2_sum(1) = x_2(1);

v_1_sum(1) = v_1(1);

v_2_sum(1) = v_2(1);

% Iterate through every frame

for j = 2:frames

x_1(2) = v_1(1);

x_2(2) = v_2(1);

v_1(2) = (-alpha.* x_1(1)) + (alpha * r_1) + (beta.* x_2(1)) – (beta.* x_1(1)) – (beta * w_1) – (beta * r_2);

v_2(2) = (-gamma.* x_2(1)) + (gamma.* x_1(1)) + (gamma * w_1) + (gamma * r_2);

% Maclaurin coefficients

for i = 2:order-1

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

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

v_1(i+1) = ((-alpha.* x_1(i)) + (beta.* x_2(i)) – (beta.* x_1(i))) / i;

v_2(i+1) = ((-gamma.*x_2(i)) + (gamma.*x_1(i))) / i;

end

% Temporary sums

x_1_sum_temp = 0;

x_2_sum_temp = 0;

v_1_sum_temp = 0;

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

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

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

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

end

% Store Maclaurin series into sum array

x_1_sum(j) = x_1_sum_temp;

x_2_sum(j) = x_2_sum_temp;

v_1_sum(j) = v_1_sum_temp;

v_2_sum(j) = v_2_sum_temp;

% Reset initial conditions

x_1(1) = x_1_sum(j);

x_2(1) = x_2_sum(j);

v_1(1) = v_1_sum(j);

v_2(1) = v_2_sum(j);

end

% Plot Results

plot(x_1_sum, x_2_sum, ‘b’);

l = legend(‘plot of x_2 sum vs x_1 sum’);

title(l,’Legend’)

title(‘Parker-Sochacki Double Spring’);

xlabel(‘x_1 mass’)

ylabel(‘x_2 mass’)

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

This simulation was also implemented in Houdini. The Parker-Sochacki method was placed into the two spheres. I then used constraints to attach the coil to the wall and spheres. I showed the first video at the top of the blog.

In these blogs, I have mentioned benefits of increasing the order of the Maclaurin series. For every video, the time step is set to 1/12 or 2/24. For this video, I lowered the order to 2 making this a similar accuracy to Verlet or Modified Euler. You will see the bobs gaining energy within a few seconds.

I could have lowered the time step but I wanted to show what happens for a higher time step at a lower order.

Here is a simulation where I gave the bob on the left an initial velocity(for the other simulations both bobs were given an initial velocity of 0) and doubled the value of the spring constant:

This has shown that at a higher order, Parker-Sochacki is handling the motion of the two bobs including their relationship to each other. Thank you for taking the time to go through this!