Greetings! In the previous blog post, I applied the Parker-Sochacki method to the Ordinary Differential Equations that make up the motion of a single spring in one dimension. I showed the math implementation on converting the system of First Order Differential Equations into a Maclaurin series for position and Maclaurin series for velocity. I showed code in Matlab and then rendered videos in Houdini.

For this post, we are following the same process applying the Parker-Sochacki method to the motion of a single spring in two dimensions. I will show how I applied the Parker-Sochacki method, code, and then rendered videos from implementing the simulation in Houdini.

Let’s start by showing the First Order Ordinary Differential equations to update position x, position y, velocity in the x direction u, and velocity in the y direction v:

I’ll break this down. k, m, x_0, y_0, b, r, and g are constants. k is the spring constant. m is the mass of the bob. x_0 and y_0 make up the position of the anchor point where the spring is suspended from. b is the damping constant. g is the constant for gravity. r is the rest length. First to make this easier, I will put constant terms together:

To apply the Parker-Sochacki method, I will introduce two auxiliary variables named a and b:

We did this so that it is possible to represent x, y, u, and v as different Maclaurin series. We will also represent a and b as a Maclaurin series:

The point of this is to represent x, y, u, v, a, and b all as polynomials that are the Maclaurin series or Taylor series centered at 0. From previous blogs I showed that if we chose to just go to the first term, this becomes a first order method. If we chose to go to the second term setting i to 2, then this becomes a second order method. If we chose to go to the 4th term setting i to 4, this becomes a fourth order method. Generally speaking, increasing the order increases the accuracy of the simulation and often leads to a larger time step than methods like Runge-Kutta 4th order, Verlet, or Forward Euler. For this simulation I chose a 9th order method.

To calculate the nth term, we need to develop a recurrence relationship between the Ordinary Differential Equations. I need to take the derivatives of auxiliary variable a and auxiliary variable b. I will also introduce a variable, c, that I don’t need the derivative of but I need to compute the Maclaurin coefficients:

We now have a system of First Order Ordinary Differential Equations. I will put these terms back into velocities u and v so you can see the whole thing together:

To compute the Maclaurin coefficients, I will need to introduce a cauchy product so that I can multiply two power series together. Also, to calculate b, I need to calculate the cube of b. I showed in the 3-body blog that you can take the power of a power series. For this implementation I will calculate the square of b and then use that to calculate the cube of b.

This means I will introduce b2 and b3.

We start with given initial conditions for x, y, u, and v. The initial condition for a is given by the formula for a above and the initial condition for b is given by the formula for b above. It is the reciprocal of a. We then use those initial conditions to create Maclaurin coefficients up to any order we desire. For this simulation, I chose 9th order.

Also, there are constants in the velocities u and v that only need to be computed when i is equal to 1. In the Maclaurin coefficients below, the i+1 for u and v only apply when i is greater than or equal to 1. This recurrence relationship comes from formulas regarding power series.

Notice that b2, b3, and c are computing the ith term as opposed to the i+1 term. This is why they do not need a Maclaurin series.

Using this formula, we now can create a Maclaurin series for x, y, u, v, a, and b:

Once the Maclaurin series is computed for every variable, they become the new initial conditions and the process repeats for as many frames as you would like.

I would also like to point out how we put everything together. If you used a numerical method like Forward Euler, Modified Euler, Verlet, or Runge-Kutta 4th order, the section that became auxiliary variable a would have to be recomputed using the square roots in the numerator and denominator with every new frame. Here everything is put together.

Next I will show code in Matlab used to implement this simulation. Because Matlab starts its index at 1, order is set to 10.

% Adam Wermus

% July 9, 2016

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

% Parker-Sochacki method

% clear everything

clear all;

clc;

% Parameters

r = 2.5; % rest length

x_0 = 0; % anchor point x

y_0 = 3; % anchor point y

g = -9.81; % gravity

k = 6.0; % spring stiffness

m = 0.5; % mass of bob

damping = 0; % damping constant

alpha = k/m;

beta = damping/m;

% Parker Sochacki setup

order = 10;% order of Maclaurin series

h = 0.04167; % time step

frames = 2500; % number of frames;

% Setup Arrays for Maclaurin coefficients

x = zeros(1, order);

y = zeros(1, order);

u = zeros(1, order);

v = zeros(1, order);

a = zeros(1, order);

b = zeros(1, order);

b2 = zeros(1, order);

b3 = zeros(1, order);

c = zeros(1, order);

% Setup Arrays that store Position and Velocity values

x_sum = zeros(1, frames);

y_sum = zeros(1, frames);

u_sum = zeros(1, frames);

v_sum = zeros(1, frames);

a_sum = zeros(1, frames);

b_sum = zeros(1, frames);

% Initial Conditions

x(1) = 0;

y(1) = -0.31667;

u(1) = 1.5;

v(1) = 1.7;

a(1) = ((x(1) – x_0).^2 + (y(1) – y_0).^2).^(1./2);

b(1) = 1./a(1);

% Initial arrays

x_sum(1) = x(1);

y_sum(1) = y(1);

u_sum(1) = u(1);

v_sum(1) = v(1);

a_sum(1) = a(1);

b_sum(1) = b(1);

% Iterate through every frame

for j = 2:frames

% 1st

c(1) = ((x(1) – x_0).*u(1)) + (y(1) – y_0).*v(1);

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

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

% Second setup

x(2) = u(1);

y(2) = v(1);

u(2) = -alpha.*x(1) + alpha.*x_0 + alpha.*r.*x(1).*b(1) – alpha.*r.*x_0.*b(1) – beta.*u(1);

v(2) = g – alpha.*y(1) + alpha.*y_0 + alpha.*r.*y(1).*b(1) – alpha.*r.*y_0*b(1) – beta.*v(1);

a(2) = c(1).*b(1);

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

% Maclaurin coefficients

for i = 2:order-1

% cauchy product for b2

cp_b2 = 0;

for k = 1:i

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

end

b2(i) = cp_b2;

% cauchy product for b3

cp_b3 = 0;

for k = 1:i

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

end

b3(i) = cp_b3;

% cauchy product for c

cp_c = 0;

for k = 1:i

cp_c = cp_c + ((x(k)-x_0).*u(i+1-k));

cp_c = cp_c + ((y(k)-y_0).*v(i+1-k));

end

c(i) = cp_c;

% Update x and y position

x(i+1) = u(i) / i;

y(i+1) = v(i) / i;

% Velocity update

cp_u = 0;

cp_v = 0;

for k = 1:i

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

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

end

cp_u = alpha.*r.*cp_u;

cp_v = alpha.*r.*cp_v;

cp_u = cp_u – (alpha.*x(i));

cp_v = cp_v – (alpha.*y(i));

cp_u = cp_u – (alpha.*r.*x_0.*b(i));

cp_v = cp_v – (alpha.*r.*y_0.*b(i));

cp_u = cp_u – (beta.*u(i));

cp_v = cp_v – (beta.*v(i));

u(i+1) = cp_u / i;

v(i+1) = cp_v / i;

% cauchy product for a and b

cp_a = 0;

cp_b = 0;

for k = 1:i

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

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

end

a(i+1) = cp_a / i;

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

end

% Temporary sums

x_sum_temp = 0;

y_sum_temp = 0;

u_sum_temp = 0;

v_sum_temp = 0;

a_sum_temp = 0;

b_sum_temp = 0;

% Compute Maclaurin series

for i = 1:order-1

x_sum_temp = x_sum_temp + (x(i)*h^(i-1));

y_sum_temp = y_sum_temp + (y(i)*h^(i-1));

u_sum_temp = u_sum_temp + (u(i)*h^(i-1));

v_sum_temp = v_sum_temp + (v(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_sum(j) = x_sum_temp;

y_sum(j) = y_sum_temp;

u_sum(j) = u_sum_temp;

v_sum(j) = v_sum_temp;

a_sum(j) = a_sum_temp;

b_sum(j) = b_sum_temp;

% Reset initial conditions

x(1) = x_sum(j);

y(1) = y_sum(j);

u(1) = u_sum(j);

v(1) = v_sum(j);

a(1) = a_sum(j);

b(1) = b_sum(j);

end

% Plot Results

plot(x_sum,y_sum, ‘k-‘);

title(‘Parker-Sochacki Single Spring 2D’)

xlabel(‘x position’);

ylabel(‘y position’);

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

Damping is set to 0 making a beautiful picture of the bob moving back and forth. It does not go to rest.

This simulation was also implemented in Houdini. In this first video, the time step is set to 1/12 and the order of the Maclaurin series is set to 9. Damping is set to zero so that the bob does not go to rest. The Parker-Sochacki implementation is in the sphere. I then used a wiresolver and static solver to constrain the coil.

For this next video, I set the damping to 0.2 so that the bob will stop moving around. I also gave it a higher initial velocity in the x direction so it begins to move to the right before constraints pull it in.

Thank you for taking the time to go through this! I hope this gives a better understanding of using the Parker-Sochacki method.