mathr / blog / #

In my previous post I tried to solve a differential equation using a port of Octave's adaptive Runge-Kutta 4/5 integration algorithm. It failed with energy explosions. Today I tried a different, much simpler, integration algorithm, and the energy explosions seem to be solved.

Velocity Verlet integration can be implemented for this problem in a few lines of code:

#include <math.h>

// compute acceleration from position
static inline void f(double *a, const double x[2])
{
  a[0] = -exp(x[1] * x[1]) * x[0];
  a[1] = -exp(x[0] * x[0]) * x[1];
}

// Velocity Verlet integration
void compute()
{
  const double h = 0.01;
  double x[2] = { 0, 1 };
  double v[2] = { 1, 1 };
  double a[2] = { 0, 0 };
  f(a, x);
  v[0] += 0.5 * h * a[0];
  v[1] += 0.5 * h * a[1];
  while (1)
  {
    x[0] += h * v[0];
    x[1] += h * v[1];
    f(a, x);
    v[0] += h * a[0];
    v[1] += h * a[1];
  }
}

You can download C99 source code for a JACK client which sonifies the chaotic coupled oscillators. I've been running it for 45mins at 48kHz sample rate, and no explosions yet, but be careful in case this is still transient behaviour...

No explosions may seem like a good thing, but what if the original differential equations are truly explosive, and the stability is a computational artifact? I still don't know the answer to that question.