Chaotic Magnetic Pendulum

See also the WebGPU version.

The plot above represents a top-down view of a chaotic magnetic pendulum. The black dots represent attracting magnets. Click and drag to move them. Each point on the plot represents the starting position of a pendulum. A pendulum is released from rest, and a bit of friction ensures it eventually comes to rest at one of the magnets. Color indicates which.

For two-dimensional position , friction , and magnets , the pendulum moves according to the equations

Pendulum height means the bottom of the pendulum is elevated slightly above the magnets so that it doesn’t experience infinite acceleration when it gets close.

For more information, see The Magnetic Pendulum or Gravity Fractals.

Implementation

The trajectory when the mouse is moved—corresponding to one trajectory—is computed on the CPU, while the colorful background field—corresponding to every trajectory—is iterated on the GPU.

In both formats, the biggest challenge is computing trajectories efficiently. The pendulum may change direction quickly and so needs a small time step, but the pendulum can take a long time to come to a rest. These competing priorities make efficient evaluation essential.

Both therefore use the adaptive RK4(5) method. The CPU variant is implemented in @rreusser/integration#ode45. Adaptive means that the timestep is modified on each step to maximize the step size while controlling the error. This is done through careful construction of the Runge-Kutta method using coefficients of the Cash-Karp method. The coefficients allow the same derivative evaluations to produce both a fifth order and so-called embedded fourth order estimate of the next state. The difference between the estimates produces an error estimate, which can subsequently be worked backwards into the timestep required to achieve a specified error threshold.

To get started, we define the derivative function

The GLSL variant of ode45 is reproduced below. It’s mostly equivalent to the CPU variant except modified for the limitations of GLSL control flow.

const float safety = 0.95;
const float maxDecrease = 0.02;
const float maxIncrease = 50.0;

vec4 ode45 (vec4 y, inout float dt) {
  // Fifth order estimate using constants for the Cash-Karp method
  vec4 k1 = deriv(y);
  vec4 k2 = deriv(y + dt * 0.2 * k1);
  vec4 k3 = deriv(y + dt * (0.075 * k1 + 0.225 * k2));
  vec4 k4 = deriv(y + dt * (0.3 * k1 - 0.9 * k2 + 1.2 * k3));
  vec4 k5 = deriv(y + dt * (
    -0.203703703703703703 * k1 +
    2.5 * k2 -
    2.592592592592592592 * k3 +
    1.296296296296296296 * k4
  ));
  vec4 k6 = deriv(y + dt * (
    0.029495804398148148 * k1 +
    0.341796875 * k2 +
    0.041594328703703703 * k3 +
    0.400345413773148148 * k4 +
    0.061767578125 * k5
  ));

  // Estimate the error using the embedded fourth order method
  vec4 tmp = dt * (
    0.004293774801587301 * k1 -
    0.018668586093857832 * k3 +
    0.034155026830808080 * k4 +
    0.019321986607142857 * k5 -
    0.039102202145680406 * k6
  );
  float err2 = dot(tmp, tmp);

  // Wasteful, but only accept the step if error is within tolerance
  bool accept = err2 <= tol2;
  if (accept) y += dt * (
    0.097883597883597883 * k1 +
    0.402576489533011272 * k3 +
    0.210437710437710437 * k4 +
    0.289102202145680406 * k6
  );

  // Either way, adjust dt according to the estimate
  dt *= clamp(
    safety * pow(tol2 / err2, accept ? 0.125 : 0.1),
    maxDecrease,
    maxIncrease
  );

  return y;
}

Note that if the error threshold is not achieved (that is, if we need to reduce the time step to keep error under control), then we simply throw out the timestep entirely. This is very wasteful! C’est la vie.

Given zero velocity and starting position xy, the final iteration loop is reproduced below. We start by packing the state into a vec4. We update the state times, overwriting the state and mutating dt on each step. The iteration is quite straightfoward, but GLSL 1.00 does not permit loop bounds to be computed at runtime, so we have to hard-code the loop extents into the shader.

vec4 y = vec4(xy, vec2(0));
float dt = 0.01;
for (int i = 0; i < ; i++) y = ode45(y, dt);

The final step is to choose a weighted average of color depending on how close it ends up to each magnet at the end of the integration loop. Has it actually come to a rest? Is it close to any of the magnets? Who knows! Or more specifically, subject to the reasonable limitations of GLSL fragment shaders, we can’t really know, so we just do our best and go with it.

I’ve chosen a rainbow color palette from Inigo Quilez and used a quasirandom sequence from Martin Roberts to distribute the colors somewhat randomly so that adjacent colors around the unit circle aren’t also adjacent around the color wheel.

That’s it. Once in place, adaptive methods are just great. It’s a pretty handy function which you’re free to lift and use. (The source is MIT Licensed.) You can view the live notebook on Observable.