Greetings! We are going to begin a series of spring simulations. I want to show how the Parker-Sochacki method handles springs.
This first blog shows how to apply the Parker-Sochacki method to the equations of motion for a single spring in one dimension. I will show the math on how I applied the Parker-Sochacki method, code in Matlab, and a rendered video in Houdini.
In a previous blog(posted in January), I solved the motion of a single 1D spring in C++ using the Parker-Sochacki method at 4th order and Runge-Kutta 4th order. Under the same time step, the same number of iterations, and the same initial conditions, Parker-Sochacki was faster but a little over 10%. I also noted that up to several decimal places, the accuracy was the same. I need to retract that statement. Unless I get the same value to all of the decimal places, the accuracy is not the same.
We will first show the Second Order Homogeneous Ordinary Differential Equation for the motion of a single spring in one dimension:
We want to solve for position and velocity. m, c, and k are constants. m represents the mass of the bob, c represents damping(friction), and k represents the stiffness of the spring.
To use the Parker-Sochacki method, we need to split this equation up into a system of First Order Ordinary Differential Equations. Then, with a given set of initial conditions for the position and velocity, we develop a recursive relationship so that we can compute any number of coefficients that will be used to create the Maclaurin series. The Parker-Sochacki method converts any system of First Order Differential Equations into a polynomial that is the Maclaurin series, or a Taylor series centered at zero. We want to represent the position, x(t), and velocity, v(t), like this:
I’ll show this step-by-step. We take the Differential Equation above and convert it into a system of First Order Ordinary Differential Equations:
Next, using properties from series solutions, we develop a recurrence relationship to compute the coefficients for the Maclaurin series:
We are given(or you choose) an x_0 and a v_0. This make it possible to always find the next coefficient. Generally speaking, the more coefficients you use for the polynomial, the more accurate the solution is(this isn’t always true as shown in previous blogs but it is true most of the time).
The t above in the Maclaurin series is the time step, which is a constant you can pick. From previous blogs, I have shown that using the Parker-Sochacki method often leads to being able to take a larger time step than other numerical methods like Runge-Kutta 2nd Order, Runge-Kutta 3rd Order, Runge-Kutta 4th Order, Forward Euler, or Verlet.
Once the Maclaurin series is evaluated, that value becomes the new x_0 for x(t) and v_0 for v(t). The process repeats for as many “frames” as you would like.
This was implemented in Matlab at 9th order. Note that if you look at the code, order will be set to 10. This is because Matlab arrays start at 1 instead of 0. When I implemented this in Houdini, I set the order constant to 9.
Here is the code implemented in Matlab. You can copy this into a matlab document, adjust the indenting, and see a plot for the position and velocity. You will see that the “bob” takes off and oscillates to rest position and a zero velocity.
% Adam Wermus
% July 6, 2016
% This code solves the motion of a single spring using the Parker-Sochacki
% Ordinary Differential Equations
% x’ = v
% v’ = (-k/m)x – (b/m)v
% Clear Content
k = 3; % spring stiffness
m = 0.5; % mass of block
b = 0.1; % damping constand(friction)
% Parker-Sochacki setup
order = 10; % order of Maclaurin series
frames = 1500; % number of frames
h = 1./24; % time step
% Setup Arrays for Maclaurin coefficients
x = zeros(1, order);
v = zeros(1, order);
% Setup Arrays that store Position and Velocity values
x_sum = zeros(1, frames);
v_sum = zeros(1, frames);
% Initial Conditions
x(1) = 0.5;
v(1) = 0.2;
x_sum(1) = x(1);
v_sum(1) = v(1);
% Iterate through every frame
for j = 2:frames
% Maclaurin coefficients
for i = 1:order-1
x(i+1) = v(i) / i;
v(i+1) = -(k/m)*(x(i)/i) – (b/m)*(v(i)/i);
% Temporary sums
x_sum_temp = 0;
v_sum_temp = 0;
% Compute Maclaurin series
for i = 1:order-1
x_sum_temp = x_sum_temp + (x(i)*h^(i-1));
v_sum_temp = v_sum_temp + (v(i)*h^(i-1));
% Store Maclaurin series into sum array
x_sum(j) = x_sum_temp;
v_sum(j) = v_sum_temp;
% Reset initial conditions
x(1) = x_sum(j);
v(1) = v_sum(j);
% Plot Results
plot(1:frames,x_sum,’b-‘, 1:frames, v_sum, ‘o’);
legend(‘position of mass’, ‘velocity of mass’);
title(‘Parker-Sochacki Single Spring’);
In Matlab, I chose a time step of 1/24.
This was also implemented in Houdini using the solver node. I have two videos where the first video has a higher initial velocity(20 m/s) than the second video(10 m/s). The time step is 2/24 or 1/12. You will see a bob move forward and then be pulled back by the spring. I applied the Parker-Sochacki method to the sphere that is the bob. Then(thank you for the help) I used constraints to attach one end of the spring to the “wall” and the other to the sphere.
Faster initial velocity:
Slower initial velocity:the sphere does not oscillate as far.
Thank you for taking the time to go through this! As a test, I set the initial velocity to 100 m/s. The “bob” where Parker-Sochacki was applied still worked properly. This means that Parker-Sochacki can handle many different initial conditions for this simulation at a time step of 1/12!