Magnetic Pendulum Basins

DEFINITION

ẍ + b·ẋ + k·x = Σₙ C·(Xₙ − x) / (|Xₙ − x|² + h²)^(5/2)
start at rest, colour by which magnet it reaches
shade by the step at which the answer stopped changing

NOTES

Hang a magnet on a string over three more arranged in a triangle, pull it aside, let go. It will end on one of the three, and asking which is a fair question with an unfair answer: the set of starting points that end on each magnet is so finely interleaved with the other two that almost everywhere on the boundary, no measurement precise enough exists. Grebogi, McDonald, Ott and Yorke gave that its name and its arithmetic in 1983 — final-state sensitivity — and showed the cost is worse than exponential: to halve your uncertainty about the outcome you must do far better than halve your uncertainty about the start, by an amount fixed by the fractal dimension of the boundary. The pendulum is the standard exhibit rather than their example, and no single first publication of the coloured picture can be identified; Peitgen, Jürgens and Saupe built a chapter on it in 1992 and it has been redrawn everywhere since. The force law here is a dipole, falling as the fourth power of distance, softened by h — the height the bob swings above the plane the magnets sit in — which is what keeps the acceleration finite directly overhead.

PROVENANCE

Origin

C. Grebogi, S. W. McDonald, E. Ott & J. A. Yorke, "Final state sensitivity: an obstruction to predictability", Physics Letters A 99(9), 1983, pp. 415–418; developed at length in S. W. McDonald, C. Grebogi, E. Ott & J. A. Yorke, "Fractal basin boundaries", Physica D 17(2), 1985, pp. 125–153. McDonald is the first author of the second and the second of the first; he is frequently dropped from both.

Exhibit

The three-magnet rendering is standard and unattributed. H.-O. Peitgen, H. Jürgens & D. Saupe, Chaos and Fractals: New Frontiers of Science, Springer, 1992, work it as their case for fractal basin boundaries; F. C. Moon, J. Cusumano & P. J. Holmes, "Evidence for homoclinic orbits as a precursor to chaos in a magnetic pendulum", Physica D 24, 1987, pp. 383–390, is the physical pendulum in the laboratory. Neither is the origin of the picture, and this plate does not claim one.

Standing

Public domain — a damped ODE and a question about its limit

Constants

k and C were scanned as a pair, because what governs the picture is their ratio and no one-at-a-time sweep finds a joint dead region: at C below 0.55 with k above 1, the integrator has not finished — 89% of starts settle at (1.20, 0.30) against 100% everywhere else — so C starts at 0.55 rather than at the 0.30 a single sweep at the default k would have allowed. The three of k, C and the view width jointly set how much of the plate is boundary rather than basin, from 14% to 74% across the declared box.

Cost

An ODE per pixel, in mandelbrot’s class with a far heavier inner loop: one whole 90-square JS pass measures 32 ms. So the JS path integrates four rows a frame into a cache and replays it — measured in a browser, no single paint exceeds 7.3 ms and the steady state once the cache is full is 0.1, which makes this the cheapest JS path in the wave rather than the most expensive. The shader is the real renderer: 0.4 ms on a grid card, and 1.2 to 3.8 ms across runs at the drawer’s full 1335×891 backing buffer, where it is the wave’s only shader that costs more than a millisecond.

Source

doi:10.1016/0375-9601(83)90945-3

TOUCHDESIGNER · GLSL

The same shader this plate runs, reframed for a GLSL TOP. Pasted bare it renders the published constants as a still frame; wire absTime.seconds into u_t on the Vectors page to animate it.

// FORMA — PL. 87 · MAGNETIC PENDULUM BASINS — Grebogi, McDonald, Ott & Yorke, 1983 · the pendulum exhibit has no single author
//   ẍ + b·ẋ + k·x = Σₙ C·(Xₙ − x) / (|Xₙ − x|² + h²)^(5/2)
//   start at rest, colour by which magnet it reaches
//   shade by the step at which the answer stopped changing
// TouchDesigner port — paste into a GLSL TOP's pixel shader. Set the
// resolution on the TOP's Common page. As pasted it renders the published
// constants as a still frame; to animate, add a uniform named u_t on the
// GLSL TOP's Vectors 1 page with the expression absTime.seconds.
// Constants are consts — edit to tweak; comments give the measured range.
// Written from the published mathematics, not adapted from any code.

#define u_res (uTDOutputInfo.res.zw)
uniform float u_t;               // absTime.seconds on the Vectors page; unset = still

const float u_phase = 0.4711;    // this plate's own grid phase, 0..1
// FORMA's FRACTALS accent as cosine-gradient coefficients
const vec3 u_pal_a = vec3(0.46, 0.1389, 0.1912);
const vec3 u_pal_b = vec3(0.5, 0.151, 0.2078);
const vec3 u_pal_c = vec3(1, 1, 1);
const vec3 u_pal_d = vec3(0, 0.05, 0.1);

const float p_b     = 0.25;        // b — damping · live 0.15 .. 0.6
const float p_k     = 0.6;         // k — restoring force · live 0.25 .. 1.2
const float p_h     = 0.4;         // h — bob height above the magnets · live 0.3 .. 0.7
const float p_mag   = 1.0;         // C — magnet strength · live 0.55 .. 2.5
const float p_span  = 3.0;         // view half-width · live 2 .. 4.5
const float p_steps = 200.0;       // integration ceiling · live 60 .. 200

/* The order's ramp — the same cosine formulation the JS kit uses, so a
   plate keeps its classification colour in either language. */
vec3 ramp(float t){
  return clamp(u_pal_a + u_pal_b * cos(6.28318530718 * (u_pal_c * t + u_pal_d)), 0.0, 1.0);
}

/* Sawtooth and triangle on this plate's phase, mirroring the JS kit. */
float cycle(float t, float period){ return fract(t / period + u_phase); }
float pingpong(float t, float period){
  float u = cycle(t, period);
  return u < 0.5 ? u * 2.0 : 2.0 - u * 2.0;
}


vec3 plate(vec2 uv){
  float TAU = 6.283185307179586;
  float ar = u_res.y / u_res.x;
  float dt = 0.06, hh = p_h * p_h;
  int cap = int(floor(p_steps + 0.5));

  float a0 = u_phase * TAU;
  vec2 m0 = vec2(cos(a0), sin(a0));
  vec2 m1 = vec2(cos(a0 + 2.0943951023931953), sin(a0 + 2.0943951023931953));
  vec2 m2 = vec2(cos(a0 + 4.1887902047863905), sin(a0 + 4.1887902047863905));

  vec2 x = vec2((uv.x - 0.5) * 2.0 * p_span, (uv.y - 0.5) * 2.0 * p_span * ar);
  vec2 v = vec2(0.0);
  int best = -1;
  float decided = 0.0;
  /* 200 is the step slider's own ceiling; the live count breaks out of it, and
     so does settling, which most pixels reach far sooner. */
  for (int s = 0; s < 200; s++){
    if (s >= cap) break;
    vec2 a = -p_k * x - p_b * v;
    vec2 d0 = m0 - x, d1 = m1 - x, d2 = m2 - x;
    float r0 = dot(d0, d0) + hh, r1 = dot(d1, d1) + hh, r2 = dot(d2, d2) + hh;
    a += d0 * (p_mag / (r0 * r0 * sqrt(r0)));
    a += d1 * (p_mag / (r1 * r1 * sqrt(r1)));
    a += d2 * (p_mag / (r2 * r2 * sqrt(r2)));
    v += a * dt;
    x += v * dt;
    /* Nearest magnet and the step it last changed on, resolved in the same
       order the JS path resolves it so the two never disagree on a tie. */
    float e0 = dot(m0 - x, m0 - x), e1 = dot(m1 - x, m1 - x), e2 = dot(m2 - x, m2 - x);
    float near = e0;
    int who = 0;
    if (e1 < near){ near = e1; who = 1; }
    if (e2 < near){ near = e2; who = 2; }
    if (who != best){ best = who; decided = float(s) + 1.0; }
    if (near < 0.05 && dot(v, v) < 0.05) break;
  }

  float age = min(1.0, decided / 120.0);
  vec3 c = ramp(0.82 + 0.15 * float(max(best, 0)));
  /* Settle time as a tone, with the sweeping band on top of it. The shader has
     no development phase because it does not need one — it computes every
     pixel every frame — so once the JS path's cache is full the two agree
     exactly: the same tone, the same band, the same settle order. */
  float lit = 0.42 + 0.45 * (1.0 - age)
            + 0.40 * max(0.0, 1.0 - abs(age - cycle(u_t, 10.0)) * 12.0);
  return c * lit;
}

out vec4 fragColor;
void main(){
  // FORMA's uv runs y-down, matching its canvas; TD's vUV runs up
  vec2 uv = vec2(vUV.s, 1.0 - vUV.t);
  fragColor = TDOutputSwizzle(vec4(plate(uv), 1.0));
}

NUKE · BLINKSCRIPT

The same shader this plate runs, transpiled to a BlinkScript kernel. Paste it into a BlinkScript node's Kernel Source and press Recompile; every constant arrives as a knob at its published value, and u_t animates with the expression frame/24. Compiled and rendered in Nuke 17.1, then compared against this plate on the page.

// FORMA — PL. 87 · MAGNETIC PENDULUM BASINS — Grebogi, McDonald, Ott & Yorke, 1983 · the pendulum exhibit has no single author
//   ẍ + b·ẋ + k·x = Σₙ C·(Xₙ − x) / (|Xₙ − x|² + h²)^(5/2)
//   start at rest, colour by which magnet it reaches
//   shade by the step at which the answer stopped changing
// Nuke port — a BlinkScript kernel. Paste into a BlinkScript node's Kernel
// Source and press Recompile. Every constant arrives as a knob at its published
// value (the comment gives the measured range); u_t is a knob too — animate it
// with the expression frame/24 or leave it at 0 for the still frame. Written
// from the published mathematics, not adapted from any code.
// Transpiled from the shader this plate runs on the page (GLSL ES 3.00):
// vec → float2/3/4, swizzles expanded, GLSL builtins Blink lacks written out
// as forma_ functions, float literals suffixed. Compiled and rendered in a
// real Nuke (17.1v1) and compared against this plate on the page: 34 of 34.
//
// plate() and its helpers are written to a single exit — the loop that runs
// once. That is not a style: Blink 17.1 drops a conditional early return from
// a called function while Vectorize is on, which is the node default, with no
// warning and no error. Written this way it paints correctly as pasted.

kernel Forma_basins : ImageComputationKernel<ePixelWise>
{
  Image<eWrite> dst;

param:
  float u_t;             // seconds; 0 is the still frame
  float p_b;     // b — damping · live 0.15 .. 0.6
  float p_k;     // k — restoring force · live 0.25 .. 1.2
  float p_h;     // h — bob height above the magnets · live 0.3 .. 0.7
  float p_mag;   // C — magnet strength · live 0.55 .. 2.5
  float p_span;  // view half-width · live 2 .. 4.5
  float p_steps; // integration ceiling · live 60 .. 200

local:
  float2 u_res;
  float u_phase;
  float3 u_pal_a, u_pal_b, u_pal_c, u_pal_d;

  void define(){
    defineParam(u_t, "u_t", 0.0f);
    defineParam(p_b, "p_b", 0.25f);
    defineParam(p_k, "p_k", 0.6f);
    defineParam(p_h, "p_h", 0.4f);
    defineParam(p_mag, "p_mag", 1.0f);
    defineParam(p_span, "p_span", 3.0f);
    defineParam(p_steps, "p_steps", 200.0f);
  }

  void init(){
    u_res = float2(float(dst.bounds.width()), float(dst.bounds.height()));
    u_phase = 0.4711f;    // this plate's own grid phase, 0..1
    // FORMA's FRACTALS accent as cosine-gradient coefficients
    u_pal_a = float3(0.46f, 0.1389f, 0.1912f);
    u_pal_b = float3(0.5f, 0.151f, 0.2078f);
    u_pal_c = float3(1.0f, 1.0f, 1.0f);
    u_pal_d = float3(0.0f, 0.05f, 0.1f);
  }

  /* GLSL builtins Blink lacks, written as templates rather than overload sets.
     Blink's operators return expression templates (Swizzle<float,N>), so a call
     passing an expression cannot resolve against an overload set on float2
     against float3 — measured in Nuke 17.1: a float2 expression is ambiguous
     between the two, while scalar-against-vector resolves. A template deduces
     the expression's own type, so the ambiguity cannot arise. */
  template <class T> T forma_fract(T v){ return v - floor(v); }
  template <class T, class S> T forma_mod(T x, S y){ return x - y * floor(x / y); }
  /* Blink's own min/max/clamp take no scalar bound against a vector, which GLSL
     does; v * 0.0f + b is that bound at the vector's own width, and collapses to
     b when v is a scalar, so one template serves both. */
  template <class T, class S> T forma_min(T a, S b){ return min(a, a * 0.0f + b); }
  template <class T, class S> T forma_max(T a, S b){ return max(a, a * 0.0f + b); }
  template <class T, class S> T forma_clamp(T v, S lo, S hi){ return clamp(v, v * 0.0f + lo, v * 0.0f + hi); }
  int forma_min(int a, int b){ return min(a, b); }
  int forma_max(int a, int b){ return max(a, b); }
  /* GLSL step(edge, x) is 1 where x >= edge; floor(sign(x - e) * 0.5 + 1) is
     that exactly, equality included, out of builtins Blink does have. */
  template <class T, class S> T forma_step(S e, T x){ return floor(sign(x - e) * 0.5f + 1.0f); }
  template <class T, class S> T forma_smoothstep(S a, S b, T x){
    T t = forma_clamp((x - a) / (b - a), 0.0f, 1.0f);
    return t * t * (3.0f - 2.0f * t);
  }
  template <class T> float forma_distance(T a, T b){ return length(a - b); }
  float forma_tanh(float x){ float e = exp(2.0f * x); return (e - 1.0f) / (e + 1.0f); }
  float forma_radians(float d){ return d * 0.01745329252f; }
  // the page's hash2 is exact uint32; Blink has int, so the shifts are made
  // logical by masking and the read-back is lifted into 0 .. 2^32
  /* A uint read back as a float. Blink has no unsigned type, so a value past
     2^31 arrives as a negative int and float() of it is negative. Measured on
     gabor, whose own generator then returned uniforms in [-0.5, 0.5) and drew
     a different picture — it compiled, it rendered, and only comparing it with

  /* The order's ramp — the same cosine formulation the JS kit uses, so a
     plate keeps its classification colour in either language. */
  float3 ramp(float t){
    return forma_clamp(u_pal_a + u_pal_b * cos(6.28318530718f * (u_pal_c * t + u_pal_d)), 0.0f, 1.0f);
  }

  /* Sawtooth and triangle on this plate's phase, mirroring the JS kit. */
  float cycle(float t, float period){ return forma_fract(t / period + u_phase); }
  float pingpong(float t, float period){
    float u = cycle(t, period);
    return u < 0.5f ? u * 2.0f : 2.0f - u * 2.0f;
  }


  float3 plate(float2 uv){
    float TAU = 6.283185307179586f;
    float ar = u_res.y / u_res.x;
    float dt = 0.06f;
    float hh = p_h * p_h;
    int cap = int(floor(p_steps + 0.5f));

    float a0 = u_phase * TAU;
    float2 m0 = float2(cos(a0), sin(a0));
    float2 m1 = float2(cos(a0 + 2.0943951023931953f), sin(a0 + 2.0943951023931953f));
    float2 m2 = float2(cos(a0 + 4.1887902047863905f), sin(a0 + 4.1887902047863905f));

    float2 x = float2((uv.x - 0.5f) * 2.0f * p_span, (uv.y - 0.5f) * 2.0f * p_span * ar);
    float2 v = float2(0.0f);
    int best = -1;
    float decided = 0.0f;
    /* 200 is the step slider's own ceiling; the live count breaks out of it, and
       so does settling, which most pixels reach far sooner. */
    for (int s = 0; s < 200; s++){
      if (s >= cap) break;
      float2 a = -p_k * x - p_b * v;
      float2 d0 = m0 - x;
      float2 d1 = m1 - x;
      float2 d2 = m2 - x;
      float r0 = dot(d0, d0) + hh;
      float r1 = dot(d1, d1) + hh;
      float r2 = dot(d2, d2) + hh;
      a += d0 * (p_mag / (r0 * r0 * sqrt(r0)));
      a += d1 * (p_mag / (r1 * r1 * sqrt(r1)));
      a += d2 * (p_mag / (r2 * r2 * sqrt(r2)));
      v += a * dt;
      x += v * dt;
      /* Nearest magnet and the step it last changed on, resolved in the same
         order the JS path resolves it so the two never disagree on a tie. */
      float e0 = dot(m0 - x, m0 - x);
      float e1 = dot(m1 - x, m1 - x);
      float e2 = dot(m2 - x, m2 - x);
      float near = e0;
      int who = 0;
      if (e1 < near){ near = e1; who = 1; }
      if (e2 < near){ near = e2; who = 2; }
      if (who != best){ best = who; decided = float(s) + 1.0f; }
      if (near < 0.05f && dot(v, v) < 0.05f) break;
    }

    float age = forma_min(1.0f, decided / 120.0f);
    float3 c = ramp(0.82f + 0.15f * float(forma_max(best, 0)));
    /* Settle time as a tone, with the sweeping band on top of it. The shader has
       no development phase because it does not need one — it computes every
       pixel every frame — so once the JS path's cache is full the two agree
       exactly: the same tone, the same band, the same settle order. */
    float lit = 0.42f + 0.45f * (1.0f - age)
              + 0.40f * forma_max(0.0f, 1.0f - fabs(age - cycle(u_t, 10.0f)) * 12.0f);
    return c * lit;
  }

  void process(int2 pos){
    // FORMA's uv runs y-down like its canvas; Nuke's rows run up
    float2 uv = float2((float(pos.x) + 0.5f) / u_res.x, 1.0f - (float(pos.y) + 0.5f) / u_res.y);
    float3 c = plate(uv);
    dst() = float4(c.x, c.y, c.z, 1.0f);
  }
};