Comparison of Integrators Using a Simple Pendulum in DynaMechs

Revision 1.0 2022-10-03

As we continue to look at the DynaMechs library, we decided to perform a quick investigation of the different integrators that it offers. DynaMechs has implementations of explicit Euler, Runge-Kutta (fourth order), and an adaptive stepsize Runge-Kutta. DynaMechs makes it easy to add more integrators, so we decided to add semi-implicit Euler, and the α-parameterized Runge-Kutta second order integrators. We then performed experiments on a simple pendulum setup in DynaMechs to see how the integrators compared; for example, the received wisdom is that explicit Euler gains energy, semi-implicit Euler conserves energy, and Runge-Kutta integrators tend to lose energy.

Pendulum Setup

To make the calculations simpler, we tried to make the DynaMechs setup as close to a simple gravity pendulum as possible. A 1mm diameter sphere of uniform density and mass 1kg was placed 1m along the X axis of a DynaMechs RevoluteJoint, the joint type chosen to constrain the sphere to move within a plane. The pendulum motion was started with an initial angle to the vertical of 0.1rad. Joint friction was set to zero and in the absence of any other resistances it was expected that the pendulum would oscillate between ±0.1 radians indefinitely.

Before simulation, the first few terms of the power series correction to the linear equation were used to obtain a prediction as to the anticipated period of the pendulum:- $$\displaystyle T=2\pi {\sqrt {\frac {L}{g}}}\left(1+{\frac {1}{16}}\theta _{0}^{2}+{\frac {11}{3072}}\theta _{0}^{4}+{\frac {173}{737280}}\theta _{0}^{6}+{\frac {22931}{1321205760}}\theta _{0}^{8}+{\frac {1319183}{951268147200}}\theta _{0}^{10} \right)$$

Which for \(L=1.0m\), \(g=9.80665ms^{-2}\) and \(\theta _{0}=0.1rad\) results in \(T=2.00766s\) (this is of course an approximation as more terms can be added to the series).

Experiments

Explicit Euler vs. Semi-Implicit Euler

The first simulation run was done with a lead-in period of 0s and a timestep of 0.00001s, with the aim of verifying the estimated period of the pendulum as calculated just above. All the integrators began simulation with the angle of the pendulum \(\theta _{0}=0.1rad\) at \(t=0s\). Looking through the dataset, we extracted the time indices when the pendulum was next at that angle (or closest to it). The outcome was that the explicit Euler was early by some margin in comparison to the other integrators, which themselves aligned very well with the estimated period:-

Integrator	Simulated pendulum period (s)
Explicit Euler	2.00320
Semi-Implicit Euler	2.00767
RK2 (Heun)	2.00767
RK2 (Midpoint)	2.00766
RK2 (Ralston)	2.00767
RK4	2.00767

To investigate explicit Euler deeper, we loaded the dataset into R to obtain a graph of the oscillations over the first 15 minutes (here the green circles show the inflection points of the pendulum movement and the horizontal lines show ±0.1 radians). It can clearly be seen that using the explicit Euler integrator the pendulum oscillations immediately begin diverging from the expected limits and the error continued to grow.

In comparison the semi-implicit Euler manages to maintain the expected ±0.1 radians:-

To continue the experiment, we reduced the simulation timestep to 0.0001s and then observed the pendulum motion for each of the integrators after simulating 1 hour, 1 day, and finally 1 month; the idea being, to allow any errors to accumulate. Looking at the explicit Euler again after only 1 hour of simulated pendulum swinging and it had diverged from any form of reality, now reaching ±0.6 radians:-

After the same lead-in time the semi-implicit Euler continued to stay constrained to ±0.1 radians:-

So it would seem that the reputations of explicit Euler to gain energy, and semi-implicit Euler to conserve energy are deserved.

Runge-Kutta

The Runge-Kutta integrators performed far better than the explicit Euler. Observing the pendulum motion after thirty days of continuous swinging the RK4 integrator still showed a pendulum period of 2.0077s, and an arc sweep of 0.2 radians (accurate to at least 9 decimal places, i.e. the pendulum was swinging between 0.100000000 rad and -0.100000000 rad). However, in an unexpected result, the various RK2 integrators we were testing (Heun, mid-point, and Ralston) all gained energy very slightly, the arc of the pendulum having very gradually increased over time. This increase was so small as to be practically negligable, however it was present. The semi-implicit Euler also seemed to gain energy, but to an even tinier degree, to the point where it might not even be a real observation:-

Integrator	Pendulum arc length at start (rad)	Pendulum arc length after thirty days (rad)
Semi-Implicit Euler	0.200000000	0.200000001
RK2 (Heun)	0.200000001	0.200000518
RK2 (Midpoint)	0.200000008	0.200000524
RK2 (Ralston)	0.200000005	0.200000519
RK4	0.200000000	0.200000000

Speed of Integrators

We would expect the calculation of the dynamics to be the major cost in an integrator, and looking at the relative durations for simulating one month of pendulum swings bears that out. The explicit and semi-implicit Euler integrators, which call the dynamics once, require half the duration of the RK2 integrators, which call the dynamics twice per timestep. The RK2 integrators are in turn half the duration of the RK4, which calls the dynamics four times per timestep. There are differences within classes of integrator, but as our computer was running other processes during the experiments it would be unreasonable to expect the durations to be practically identical.

Integrator	Time take to simulate 1 month (2,592,020s) of pendulum swings
Explicit Euler	490.189s (8m 10.2s)
Semi-Implicit Euler	468.467s (7m 48.5s)
RK2 (Heun)	925.802s (15m 25.8s)
RK2 (Midpoint)	929.525s (15m 29.5s)
RK2 (Ralston)	929.627s (15m 29.6s)
RK4	1862.450s (31m 2.5s)

Conclusions

This was just a quick look at some integrators available with the DynaMechs library, with experiments on the simplest possible mechanical setup, a simple gravity pendulum. We were expecting the explicit Euler to not perform particularly well, but the rate at which it deteriorated was a surprise — even using a very small time step it would seem to be a poor choice. The semi-implicit Euler performs far better for that kind of oscillatory motion, without requiring more calls to the dynamics.

Of course, the choice of integrator very much depends on the nature of the system and the final deliverable. For the kind of simple linkage systems and short timescales we are experimenting with (robot arms etc.) most integrators are probably fine. The idea that RK4 loses energy is probably from applications like orbital dynamics where planetary systems are simulated over very long timescales and to very high degrees of precision. The observation that the RK2 integrators might be gaining energy is a surprise and might require a further look when we have the time.

Source Code

/**** Integrator pendulum test *****

g++ -c pendulum.cpp -I/usr/local/include/AS-DynaMechs/ -DDM_DOUBLE_PRECISION -D__sgi -pedantic -Wpointer-arith -O2 -s -D_BOOL -Wuninitialized -c -fPIC -D_POSIX_C_SOURCE=199506L
g++ -o pendulum pendulum.o -lASdm -ldmu -fPIC -lm -ldl -lrt
*/

#include <iostream>
#include <fstream>

#include <dmTime.h>
#include <dm.h> 

#include <dmSystem.hpp>
#include <dmArticulation.hpp>
#include <dmLink.hpp>
#include <dmRevoluteLink.hpp>

#include <dmEnvironment.hpp>

#include <dmIntegrator.hpp>
#include <dmIntegSemiImplicitEuler.hpp>	// new integrator
#include <dmIntegEuler.hpp>
#include <dmIntegRK2.hpp>					// new integrator
#include <dmIntegRK4.hpp>

// compute a more accurate period for a simple pendulum using the first 
// few terms of the power series correction to the linear equation
Float simple_pendulum_power_series( Float theta, Float L )
{
	Float g = 9.80665;
	Float c[6] = { 1.0, 1.0/16.0, 11.0/3072.0, 173.0/737280.0, 22931.0/1321205760.0, 1319183.0/951268147200.0 };

	Float initial = 2*M_PI*sqrt(L/g);
	Float series = 1.0;
	for ( unsigned int i = 1; i <= 5; i++ ) {
		series += c[i] * pow(theta, i*2);
	}
	
	return initial*series;
}

int main(int argc, char** argv)
{

	CartesianVector gravity = { 0.0, 0.0, (Float)-9.80665 };
	dmEnvironment *environment = new dmEnvironment();
	dmEnvironment::setEnvironment(environment);
	environment->setGravity(gravity);

	dmArticulation *Pendulum;
	dmRevoluteLink *r1;

    Pendulum = new dmArticulation;	

	CartesianVector pos = {0.0, 0.0, 10.0};
    Quaternion quat = {0.0, 0.0, 0.0, 1.0};
   
    Pendulum->setRefSystem(quat, pos);

	/**** R1 - pendulum body ******/
	r1 = new dmRevoluteLink();	

	Float r1_mass = 1.0;
	CartesianVector r1_CoG = { 1.0, 0.0, 0.0 };
	Float r1_Ixx = 0.0000001;
	Float r1_Iyy = 1.0000001;
	Float r1_Izz = 1.0000001;
	CartesianTensor r1_I_bar = {	{r1_Ixx,   0.0,		0.0},
			                        {0.0,      r1_Iyy,	0.0},
                            		{0.0,      0.0,		r1_Izz}	};

    r1->setInertiaParameters(r1_mass, r1_I_bar, r1_CoG);    
    r1->setMDHParameters(0.0, M_PI/2.0, 0.0, 0.0);
    r1->setJointFriction( 0.0 );
	
    Pendulum->addLink(r1, NULL);

	/*********/

	/* Setup integrators */
	enum Integrators { EXPLICIT_EULER, SEMI_IMPLICIT_EULER, HEUN_RK2, MIDPOINT_RK2, RALSTON_RK2, RK4, NUM_INTEGRATORS };

	dmIntegrator *integrators[NUM_INTEGRATORS];
	integrators[EXPLICIT_EULER] 		= new dmIntegEuler();
    integrators[SEMI_IMPLICIT_EULER] 	= new dmIntegSemiImplicitEuler();
	integrators[HEUN_RK2] 				= new dmIntegRK2(RK2_HEUN);
	integrators[MIDPOINT_RK2] 			= new dmIntegRK2(RK2_MIDPOINT);
	integrators[RALSTON_RK2] 			= new dmIntegRK2(RK2_RALSTON);
	integrators[RK4] 					= new dmIntegRK4();

	// for output table columns
	string integrator_names[NUM_INTEGRATORS] = { "ExplicitEuler", "SemiImplicitEuler", "RK2Heun", "RK2Midpoint", "RK2Ralston", "RK4" };

	Float idt = 0.00001;						// simulation timestep
	const Float record_duration = 20.0;   	// seconds for which to record joint states (after lead-in)
	const unsigned int step = 1;          	// sampling rate for data output
	const unsigned int max_steps = (int)(record_duration / idt);
	const unsigned int max_record = max_steps / step;
	const Float lead_in = 0;				// seconds to run before recording states
	const unsigned int lead_in_steps = (int)(lead_in / idt);
	Float *States[NUM_INTEGRATORS];
    Float joint_pos_s = 0.1+(-M_PI/2.0);
	Float joint_vel_s = 0.0;
    Float joint_pos, joint_vel;

	Float period = simple_pendulum_power_series( 0.1, 1.0 );
	cout << "Expected period of pendulum:\t" << period << endl;
	
	for ( unsigned int i = 0; i < NUM_INTEGRATORS; i++ ) {

		cout << "Running: " << integrator_names[i] << "... ";

        // reset pendulum, set up integrator
        r1->setState(&joint_pos_s, &joint_vel_s);
        integrators[i]->addSystem(Pendulum);
		States[i] = new Float[max_record];
    
		// time the simulation
		dmTimespec start_tv, end_tv;
		dmGetSysTime(&start_tv);

		// simulate for the lead-in time, errors accumulate
    	for ( unsigned int j = 0; j < lead_in_steps; j++ ) 
		{
            integrators[i]->simulate( idt );
		}

		// record joint state for a time after lead-in
		// we will just record/sample all the states and 
		// post-process the data later with R
		unsigned int k = 0;
	
		for ( unsigned int j = 0; j < max_steps; j++ ) {

			// if necessary sample it so we don't use too much memory
			if ( (j % step == 0) && (k < max_record) ) {

                r1->getState(&joint_pos, &joint_vel);
				States[i][k] = joint_pos+(M_PI/2.0);	// adjust the angle so that the pendulum oscillates around 0 instead of pi/2
				k++;
			}

            integrators[i]->simulate( idt );
		}

		dmGetSysTime(&end_tv);
		double elapsed_time = ((double) end_tv.tv_sec - start_tv.tv_sec) + (1.0e-9*((double) end_tv.tv_nsec - start_tv.tv_nsec));

		cout << "Done. Simulated " << record_duration+lead_in << " seconds in " << elapsed_time << " seconds." << endl;	
	}

	/* write data to file for analysis */
	fstream out_file;
	out_file.open("pendulum_data.dat", ios::out);
	if ( !out_file ) {
		cerr << "Output file could not be created" << endl;
		return 1 ;
	}
	
	// write header
	out_file << "Time\tLimit1\tLimit2";
	for ( unsigned int i = 0; i < NUM_INTEGRATORS; i++ ) {
		out_file << "\t" << integrator_names[i];
	}
	out_file << endl;

	// write data
	Float timestep = 0.0;
	out_file << fixed << setprecision(12);
	for ( unsigned int j = 0; j < max_record; j++ ) {
		
		out_file << timestep << "\t0.1\t-0.1";
		for ( unsigned int i = 0; i < NUM_INTEGRATORS; i++ ) {
			out_file << "\t" << States[i][j];
		}
		out_file << endl;
		
		timestep += idt * step;
	}
	out_file.close();

/* Example of producing a graph of the data using R:-
> # load the data
> pend_data <- read.table("pendulum_data.dat", header=TRUE, sep="\t", quote="", dec=".");
> # differentiate to find inflection points (where the pendulum oscillation turns over)
> SIEulerDiff <- c(FALSE, diff(diff(pend_data$SemiImplicitEuler)>0)!=0)
> # main plot
> plot(pend_data$Time[SIEulerDiff], pend_data$SemiImplicitEuler[SIEulerDiff], type="l", col="3", main="Semi-Implicit Euler", xlab="Time (s) + 0s lead-in", ylab="Joint angle (rad)")
> # add lines for expected oscillation limits
> lines(pend_data$Time, pend_data$Limit1, type="l", col="1")
> lines(pend_data$Time, pend_data$Limit2, type="l", col="1")
> # save the plot to a file
> savePlot(filename="pendulum.png", type="png")
*/

	return 0;
}