Lyapunov Fractal

DEFINITION

xₙ₊₁ = rₙ·xₙ(1 − xₙ),  rₙ ∈ {a, b} by a repeating word
λ = lim (1/N)·Σ ln|rₙ(1 − 2xₙ)|
λ < 0 stable · λ > 0 chaotic

NOTES

Drive the logistic map with two growth rates instead of one, alternating between them on a fixed repeating word, and plot the Lyapunov exponent over the (a, b) plane. Negative means nearby starts converge and the population settles; positive means they separate and it never does. Markus and Hess found the boundary between those regions is not a curve but an intricate structure of interlocking arms, and the shape depends entirely on the word — swapping AB for AABB rebuilds the picture.

PROVENANCE

Origin

M. Markus & B. Hess, "Lyapunov Exponents of the Logistic Map with Periodic Forcing", Computers & Graphics 13(4), 1989

Standing

Public domain — an exponent of a published map

Cost

Per-pixel iteration, in mandelbrot’s class. Iteration count is locked, and the shader is why this renders sharp.

Source

doi:10.1016/0097-8493(89)90019-8

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. 65 · LYAPUNOV FRACTAL — Mario Markus & Benno Hess, 1989
//   xₙ₊₁ = rₙ·xₙ(1 − xₙ),  rₙ ∈ {a, b} by a repeating word
//   λ = lim (1/N)·Σ ln|rₙ(1 − 2xₙ)|
//   λ < 0 stable · λ > 0 chaotic
// 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.9983;    // 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_seq      = 0.0;         // forcing word · live 0 .. 3
const float p_iter     = 90.0;        // iterations · live 40 .. 120
const float p_contrast = 2.6;         // contrast · live 1 .. 6
const float p_lo       = 2.5;         // window — low rate · live 2 .. 3.8
const float p_span     = 1.5;         // window — span · live 0.2 .. 2

/* 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){
  int s = int(floor(p_seq + 0.5));
  int bits = 2, len = 2;                        // AB
  if (s == 1){ bits = 12; len = 4; }            // AABB
  else if (s == 2){ bits = 4; len = 3; }        // AAB
  else if (s == 3){ bits = 14; len = 4; }       // ABBB
  int n = int(floor(p_iter + 0.5));
  float lo = p_lo, hi = min(4.0, p_lo + p_span);
  float a = lo + (hi - lo) * uv.x;
  float b = lo + (hi - lo) * (1.0 - uv.y);
  float x = 0.5, sum = 0.0;
  for (int i = 0; i < 120; i++){                // 120 is the slider's own ceiling
    if (i >= n) break;
    float r = ((bits >> (i % len)) & 1) == 1 ? b : a;
    x = r * x * (1.0 - x);
  }
  for (int i = 0; i < 120; i++){
    if (i >= n) break;
    float r = ((bits >> ((i + n) % len)) & 1) == 1 ? b : a;
    x = r * x * (1.0 - x);
    float d = abs(r * (1.0 - 2.0 * x));
    sum += log(max(d, 1e-12));
  }
  float lam = sum / float(n);
  if (lam > 0.0) return ramp(0.52);
  return ramp(0.86 + 0.22 * min(1.0, -lam * p_contrast));
}

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. 65 · LYAPUNOV FRACTAL — Mario Markus & Benno Hess, 1989
//   xₙ₊₁ = rₙ·xₙ(1 − xₙ),  rₙ ∈ {a, b} by a repeating word
//   λ = lim (1/N)·Σ ln|rₙ(1 − 2xₙ)|
//   λ < 0 stable · λ > 0 chaotic
// 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_lyapunov : ImageComputationKernel<ePixelWise>
{
  Image<eWrite> dst;

param:
  float u_t;             // seconds; 0 is the still frame
  float p_seq;      // forcing word · live 0 .. 3
  float p_iter;     // iterations · live 40 .. 120
  float p_contrast; // contrast · live 1 .. 6
  float p_lo;       // window — low rate · live 2 .. 3.8
  float p_span;     // window — span · live 0.2 .. 2

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_seq, "p_seq", 0.0f);
    defineParam(p_iter, "p_iter", 90.0f);
    defineParam(p_contrast, "p_contrast", 2.6f);
    defineParam(p_lo, "p_lo", 2.5f);
    defineParam(p_span, "p_span", 1.5f);
  }

  void init(){
    u_res = float2(float(dst.bounds.width()), float(dst.bounds.height()));
    u_phase = 0.9983f;    // 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){
    float3 forma_r = float3(0.0f, 0.0f, 0.0f);
    for (int forma_once = 0; forma_once < 1; forma_once++){
      int s = int(floor(p_seq + 0.5f));
      int bits = 2, len = 2;                        // AB
      if (s == 1){ bits = 12; len = 4; }            // AABB
      else if (s == 2){ bits = 4; len = 3; }        // AAB
      else if (s == 3){ bits = 14; len = 4; }       // ABBB
      int n = int(floor(p_iter + 0.5f));
      float lo = p_lo;
      float hi = forma_min(4.0f, p_lo + p_span);
      float a = lo + (hi - lo) * uv.x;
      float b = lo + (hi - lo) * (1.0f - uv.y);
      float x = 0.5f;
      float sum = 0.0f;
      for (int i = 0; i < 120; i++){                // 120 is the slider's own ceiling
        if (i >= n) break;
        float r = ((bits >> (i % len)) & 1) == 1 ? b : a;
        x = r * x * (1.0f - x);
      }
      for (int i = 0; i < 120; i++){
        if (i >= n) break;
        float r = ((bits >> ((i + n) % len)) & 1) == 1 ? b : a;
        x = r * x * (1.0f - x);
        float d = fabs(r * (1.0f - 2.0f * x));
        sum += log(forma_max(d, 1e-12));
      }
      float lam = sum / float(n);
      if (lam > 0.0f) { forma_r = ramp(0.52f); break; }
      { forma_r = ramp(0.86f + 0.22f * forma_min(1.0f, -lam * p_contrast)); break; }
    }
    return forma_r;
  }

  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);
  }
};