Fractal Geometry — Mandelbrot set, Hausdorff dimension and iterated function systems

1. Self-similarity and fractal dimension

A fractal is a set (or geometric object) that exhibits self-similarity: zooming in reveals structure that resembles the whole. The key mathematical signature is a non-integer Hausdorff dimension — a fractal lives "between" classical dimensions.

For a line: N=2 copies at r=1/2, so D = log(2)/log(2) = 1. ✔
For a square: N=4, r=1/2, D = log(4)/log(2) = 2. ✔
For the Koch snowflake: N=4, r=1/3, D = log(4)/log(3) ≈ 1.262. A curve with dimension >1.

2. The Mandelbrot set and escape-time algorithm

The Mandelbrot set M is the set of complex numbers c for which the iteration zₙ₊₁ = zₙ² + c (starting from z₀ = 0) remains bounded — that is, |zₙ| never exceeds the escape radius (typically 2).

The "escape time" — the number of iterations before |z| > 2 — is used to map each pixel to a colour. Points that never escape (those in M) are traditionally black; the surrounding boundary creates the infinitely detailed fringe.

// Mandelbrot escape-time — CPU reference implementation
function mandelbrot(cx, cy, maxIter = 256) {
  let x = 0, y = 0;
  for (let i = 0; i < maxIter; i++) {
    const x2 = x * x, y2 = y * y;
    if (x2 + y2 > 4) return i;   // escaped — return iteration count
    y = 2 * x * y + cy;
    x = x2 - y2 + cx;
  }
  return maxIter;  // did not escape — inside the set
}

// Render to Canvas — maps viewport [x0,x1]×[y0,y1] to canvas pixels
function render(ctx, w, h, x0, x1, y0, y1) {
  const img = ctx.createImageData(w, h);
  for (let py = 0; py < h; py++) {
    for (let px = 0; px < w; px++) {
      const cx = x0 + (px / w) * (x1 - x0);
      const cy = y0 + (py / h) * (y1 - y0);
      const n  = mandelbrot(cx, cy);
      const t  = n / 256;  // 0..1
      const i4 = (py * w + px) * 4;
      // Simple blue-to-gold palette
      img.data[i4]   = Math.sin(t * 3.14159        ) * 255 | 0;
      img.data[i4+1] = Math.sin(t * 3.14159 * 2) * 255 | 0;
      img.data[i4+2] = t < 0.5 ? t * 2 * 255 | 0 : 255;
      img.data[i4+3] = 255;
    }
  }
  ctx.putImageData(img, 0, 0);
}

3. Smooth colouring — the continuous iteration count

Integer escape counts produce visible colour bands ("lava lamp" banding). The continuous (smooth) iteration count removes this artefact by using the final |z| value to interpolate between integer bands:

// Smooth iteration count (CPU)
function mandelbrotSmooth(cx, cy, maxIter = 512) {
  let x = 0, y = 0;
  for (let i = 0; i < maxIter; i++) {
    const x2 = x*x, y2 = y*y;
    if (x2 + y2 > 256) {
      // Smooth escape (escape radius 16, bailout at 256 for better smoothing)
      return i - Math.log2(Math.log2(Math.sqrt(x2 + y2)));
    }
    y = 2*x*y + cy;
    x = x2 - y2 + cx;
  }
  return maxIter;
}

// Smooth palette lookups (HSL cycling)
function palette(t) {
  const h = (t * 360 * 0.15) % 360;
  return `hsl(${h}, 80%, 55%)`;
}

4. Julia sets and parameter space

A Julia set Jc is computed with the same iteration zₙ₊₁ = zₙ² + c but this time c is fixed and the starting value z₀ varies over every pixel. The filled Julia set Kc is the set of z₀ that do not escape.

The Mandelbrot set is the "map" of all Julia sets: each point c in M corresponds to a connected Julia set; each point outside M corresponds to a disconnected (Cantor-like dust) Julia set. Zooming into the Mandelbrot boundary at a point c reveals a local self-similar copy of the corresponding Julia set Jc.

Computing a Julia set is identical to the Mandelbrot loop but with z₀ = pixel coordinates and c as a uniform parameter — trivial to add to any Mandelbrot shader with one boolean toggle.

5. Iterated Function Systems (IFS)

An IFS is a finite collection of contraction mappings {f₁, f₂, …, fₙ} on ℝ² (or ℝ³). By the Banach fixed-point theorem, the unique compact set that is invariant under all mappings simultaneously — the attractor — is a fractal.

The random iteration algorithm (Chaos Game) constructs the attractor by starting from any point and repeatedly applying a randomly chosen mapping. The orbit converges to the attractor after a few dozen warm-up steps.

// Barnsley fern IFS (4 affine maps, chosen probabilistically)
const fern = [
  // [a, b, c, d, e, f, prob]  — transforms [x,y] → [ax+by+e, cx+dy+f]
  [  0.00,  0.00,  0.00,  0.16,  0,     0,     0.01],  // stem
  [  0.85,  0.04, -0.04,  0.85,  0,     1.60,  0.85],  // large leaflets
  [  0.20, -0.26,  0.23,  0.22,  0,     1.60,  0.07],  // small left leaflets
  [ -0.15,  0.28,  0.26,  0.24,  0,     0.44,  0.07],  // small right leaflets
];

function chaosGame(points = 200_000) {
  let x = 0, y = 0;
  const out = [];
  for (let i = 0; i < points; i++) {
    const r = Math.random();
    let cumP = 0, fi;
    for (fi of fern) { cumP += fi[6]; if (r < cumP) break; }
    const [a,b,c,d,e,f] = fi;
    const nx = a*x + b*y + e;
    const ny = c*x + d*y + f;
    x = nx; y = ny;
    out.push(x, y);
  }
  return out;
}

Sierpiński triangle via Chaos Game

// Three vertices of an equilateral triangle
const verts = [[0.5,0], [0,1], [1,1]];
let [x, y] = [0.5, 0.5];
for (let i = 0; i < 100_000; i++) {
  const [vx, vy] = verts[Math.floor(Math.random() * 3)];
  x = (x + vx) / 2;   // move halfway to chosen vertex
  y = (y + vy) / 2;
  plotPixel(x, y);       // the Sierpiński triangle emerges
}

6. Newton fractals and basin attraction

Apply Newton's root-finding method to a polynomial p(z) in the complex plane. For each pixel c = z₀, iterate until z converges to a root:

The boundaries between basins of attraction are fractals — the Julia set of the Newton map. For polynomials of degree ≥ 3, the boundary is densely interleaved and fractal between every pair of roots.

// Newton fractal for z³ - 1 = 0
function newton(cx, cy, maxIter = 64) {
  let x = cx, y = cy;
  for (let i = 0; i < maxIter; i++) {
    const r2 = x*x + y*y, r6 = r2*r2*r2;
    if (r6 < 1e-15) break;       // at z=0 → derivative 0 (rare)
    // p(z) = z³ - 1,  p'(z) = 3z²
    // z' = z - (z³-1)/(3z²) = (2z³+1)/(3z²)
    const denom = 3 * r2 * r2;
    // numerator: 2z³+1 — compute z³ first
    const z3x = x*(x*x - 3*y*y), z3y = y*(3*x*x - y*y);
    const nx = (2*z3x + 1) / (3*r2);
    const ny = (2*z3y)      / (3*r2);
    // compute (2z³+1) / (3z²) via complex division
    x = (nx * x + ny * y) / r2;  // simplified: Re[(2z³+1)/(3z²)]
    y = (ny * x - nx * y) / r2;
  }
  // Identify closest root and return (root index, iteration count)
  const roots = [[1,0],[-0.5,0.866],[-0.5,-0.866]];
  let minD = Infinity, rootIdx = 0;
  for (let r = 0; r < 3; r++) {
    const d = (x-roots[r][0])**2 + (y-roots[r][1])**2;
    if (d < minD) { minD = d; rootIdx = r; }
  }
  return rootIdx;
}

7. GPU rendering in WebGL2 — deep zoom at 60 fps

CPU rendering at 800×600 × 512 iterations is about 250 million operations per frame — far too slow for interactive zoom. A fragment shader computes the escape-time iteration for every pixel in parallel on the GPU.

// Mandelbrot GLSL fragment shader (WebGL2)
#version 300 es
precision highp float;
uniform vec4 uView;      // (x0, y0, width, height) in complex plane
uniform int  uMaxIter;
in vec2  vUv;
out vec4 fragColor;

// Smooth iteration count + palette
vec3 palette(float t) {
  vec3 a = vec3(0.5), b = vec3(0.5);
  vec3 c = vec3(1.0);
  vec3 d = vec3(0.263, 0.416, 0.557);
  return a + b * cos(6.28318 * (c * t + d));  // Iq cosine palette
}

void main() {
  float cx = uView.x + vUv.x * uView.z;
  float cy = uView.y + vUv.y * uView.w;
  float x = 0.0, y = 0.0;
  int   n = 0;

  for (int i = 0; i < 512; i++) {
    if (i >= uMaxIter) break;
    float x2 = x*x, y2 = y*y;
    if (x2 + y2 > 256.0) { n = i; break; }
    y = 2.0*x*y + cy;
    x = x2 - y2 + cx;
    n++;
  }

  if (n == uMaxIter) { fragColor = vec4(0.0, 0.0, 0.0, 1.0); return; }

  // Smooth colouring
  float logZn = log(x*x + y*y) * 0.5;
  float nu    = log(logZn / log(2.0)) / log(2.0);
  float t     = (float(n) + 1.0 - nu) / float(uMaxIter);
  fragColor = vec4(palette(t), 1.0);
}

8. Fractals in nature, science and games

• Terrain generation: fractional Brownian motion (fBm) — layers of Perlin noise at different frequencies — produces landscapes with fractal character. The Hurst exponent H controls roughness; H=0.5 gives Brownian noise, H→1 smooth.
• Coastlines and river networks: Richardson (1961) showed that measured coastline length grows without bound as ruler length shrinks. River drainage networks show fractal branching with D ≈ 1.2.
• Lung and vascular trees: bronchial trees and blood-vessel networks are self-similar fractal trees. The fractal branching maximises surface area (gas exchange / nutrient delivery) within a bounded volume.
• Antenna design: fractional IFS antennas (Sierpiński triangle, Koch curve) resonate at multiple self-similar frequencies — compact wideband antennas used in mobile phones.
• Image compression (IFS coding): Michael Barnsley's fractal image compression finds the IFS whose attractor best approximates a given image. Decompression is free (iterate the IFS); compression is expensive (search problem).
• Game procedural generation: Midpoint displacement and Diamond-Square algorithms generate height maps for terrains using recursive self-similar subdivision. The Mandelbrot set is used for procedural biome shading in several indie games.
• Financial time series: Mandelbrot (1963) modelled cotton prices with stable Lévy distributions — the price charts exhibit statistical self-similarity across time scales (fractal market hypothesis).