This blog shows how to apply the Parker-Sochacki method to the motion of a spring subject to damping. I show code on how to implement this using Parker-Sochacki and then code on how to implement this using Runge-Kutta 4th Order in C++. Then I discuss how the two methods compare in terms of speed.

We discuss the following Ordinary Differential Equation:

where m, c, and k are constants. m is the mass, c is a damping constant, and k is a spring constant. This Second Order Ordinary Differential Equation represents the motion.

I wrote in a pdf how to apply Parker-Sochacki. This ODE needs to be split into two First Order Ordinary Differential Equations and then develop a recurrence relationship:

x represents position and v represents velocity. x(0) and v(0) are given so that the next coefficients can be computed.

In previous posts, I showed code in Matlab that used arrays to store Maclaurin coefficients so that you could adjust the order of Parker-Sochacki. For this example, I am strictly focusing on Parker-Sochacki 4th Order. I hard coded it. This means that the accuracy between Runge-Kutta 4th Order and Parker-Sochacki 4th Order should be the same under the same time step. It is. First is the Runge-Kutta implementation. This program solves RK4 and then writes the data out to a file. I also have a timer to see how long this simulation took.

C++ Code: Runge-Kutta 4th Order for a mass spring model:

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

/* Adam Wermus

April 7, 2015

Initial Value Problems with RK4

*/

#include <iostream>

#include <fstream>

#include <vector>

#include <string>

#include <time.h>

using namespace std;

// declare functions

// dx function

double dx(double v)

{

double dx;

dx = v;

return dx;

}

// dv function

double dv(double x, double v)

{

double dv;

double k = 5; // spring constant

double c = 4; // damping coefficient

double m = 5; // mass

dv = (-c*v – k*x) / m;

return dv;

}

//Runge Kutta 4th Order Function

void rk4(double x, double v, double h, ofstream& outfile)

{

double k1, k2, k3, k4, l1, l2, l3, l4;

for (double i = 0; i <= 100; i = i + h)

{

// k and l values

k1 = dx(v);

l1 = dv(x, v);

k2 = dx(v + l1*h / 2);

l2 = dv(x + k1*h / 2, v + l1*h / 2);

k3 = dx(v + l2*h / 2);

l3 = dv(x + k2*h / 2, v + l2*h / 2);

k4 = dx(v + l3*h);

l4 = dv(x + k3*h, v + l3*h);

// Outfile

outfile << “\t” << i << “\t” << x << “\t” << v << endl;

// x value

x = x + h / 6.0 * (k1 + 2 * k2 + 2 * k3 + k4);

// v value

v = v + h / 6.0 * (l1 + 2 * l2 + 2 * l3 + l4);

}

}

int main()

{

clock_t t1, t2; // timer begin

t1 = clock();

// declare variables

double h = 0.05;

string filename;

int i = 0;

// initial conditions

double x = 0.3;

double v = 0;

cout << “Enter the name of the output file” << endl;

cin >> filename;

ofstream outfile;

outfile.open(filename.c_str());

rk4(x, v, h, outfile);

outfile.close();

t2 = clock();

float diff = ((float)t2 – (float)t1);

float seconds = diff / CLOCKS_PER_SEC;

cout << seconds << endl;

system(“PAUSE”);

return 0;

}

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

Next is the Parker-Sochacki implementation:

C++ Code to Solve Parker-Sochacki 4th Order:

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

/* Adam Wermus

January 13, 2016

This program solves the mass-spring system ODE using Parker-Sochacki at order 4

x'(t) = v

v'(t) = (-cv – kx) / m

where

c = 4 // Damping Constant

k = 5 // Spring Constant

m = 5 // mass

*/

#include <iostream>

#include <fstream>

#include <vector>

#include <string>

#include <time.h>

using namespace std;

int main() {

clock_t t1, t2; // timer begin

t1 = clock();

// Setup constant values

const float c = 4;

const float k = 5;

const float m = 5;

const float inv_mass = 1 / m;

// Setup x_sum and v_sum

float x_sum = 0;

float v_sum = 0;

// Setup initial conditions

float x_0 = 0.3;

float v_0 = 0;

// Setup Maclaurin Coefficients

float x_1 = 0;

float v_1 = 0;

float x_2 = 0;

float v_2 = 0;

float x_3 = 0;

float v_3 = 0;

float x_4 = 0;

float v_4 = 0;

// Setup time step

const float dt = 0.05;

float dt2 = dt * dt;

float dt3 = dt2 * dt;

float dt4 = dt3 * dt;

// filename

string filename;

cout << “enter the name of the output file” << endl;

cin >> filename;

ofstream outfile;

outfile.open(filename.c_str());

// print first part

x_sum = x_0;

v_sum = v_0;

float check = 0;

for (int i = 0; i <= 1980; i++)

{

// Maclaurin coefficients

x_1 = v_0;

v_1 = (-c * v_0 – k * x_0) * inv_mass;

x_2 = v_1 / 2;

v_2 = (-c * v_1 – k * x_1) * inv_mass * 0.5;

x_3 = v_2 / 3;

v_3 = (-c * v_2 – k * x_2) * inv_mass * 0.3333333333;

x_4 = v_3 / 4;

v_4 = (-c * v_3 – k * x_3) * inv_mass * 0.25;

// x_sum and v_sum

x_sum = x_0 + (x_1 * dt) + (x_2 * dt2) + (x_3 * dt3) + (x_4 * dt4);

v_sum = v_0 + (v_1 * dt) + (v_2 * dt2) + (v_3 * dt3) + (v_4 * dt4);

// Outfile

outfile << “\t” << i << “\t” << x_sum << “\t” << v_sum << endl;

// update check

check += 0.05;

// reset initial conditions

x_0 = x_sum;

v_0 = v_sum;

// reset sums

x_sum = 0;

v_sum = 0;

}

outfile.close();

t2 = clock();

float diff = ((float)t2 – (float)t1);

float seconds = diff / CLOCKS_PER_SEC;

cout << seconds << endl;

system(“pause”);

return 0;

}

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

The first thing I want to show is this. For this example, Parker-Sochacki at 4th Order achieved the same accuracy as RK4 for Position and Velocity. This was expected but I wanted to confirm this seeing it in action. After 200 iterations at a time step of 0.05. The Position and Velocity data were as follows:

Runge-Kutta 4th Order Position: -0.00461

Runge-Kutta 4th Order Velocity: -0.00184

Parker-Sochacki 4th Order Position: -0.00461

Parker-Sochacki 4th Order Velocity: -0.00184

The next test was testing how long both programs took to run. Granted, these solutions might be done in such a way that I overlooked why one method was faster. But the speed was very close. I set both to more iterations than what is seen in the code because the timer sometimes registered both as running for zero seconds. I ran both programs ten times and took the mean value. Runge-Kutta 4th Order performed at a mean value of 0.006 seconds and Parker-Sochacki 4th Order performed at a mean value of 0.0053 seconds making it slightly faster. However, these values are very close and could have been different if implemented again. Also Parker-Sochacki had the highest runtime of 0.01 for one case where Runge-Kutta’s highest runtime was 0.009 seconds.

Another note on this. This does not mean Parker-Sochacki at 4th Order is always faster than Runge-Kutta 4th Order. If this system of Ordinary Differential Equations required using cauchy products, that would have most likely slowed down the simulation.

There is one serious advantage to using Parker-Sochacki that I consistently use in previous blogs but did not use for this test. The order of Parker-Sochacki can change. Parker-Sochacki could have changed to a 5th order or 6th order or 10th order Maclaurin series improving its accuracy over Runge-Kutta 4th Order. This would slow down the simulation but would open up the possibility of a larger time step. That is why Parker-Sochacki is so powerful.

If this was being implemented into a film or game, I would have raised the order of Parker-Sochacki playing with the time step to see what is the optimal way to use Parker-Sochacki for that particular simulation.