The Newton Fractal

If we apply Newton's root-finding method to a polynomial in the complex plane and color each starting point by which root it converges to and how quickly it arrives, the boundary between the root basins unfolds into an intricate fractal known as the Newton fractal.

World overview

For the cubic polynomial p(z) = z³ − 1, the three cube roots of unity sit on the unit circle at 120° intervals: r0 = 1 on the positive real axis and r1, r2 = −1/2 ± i√3/2 in the left half-plane.

In 1879, Arthur Cayley asked how Newton's method partitions the complex plane among the three cube roots of unity. While each root is surrounded by a broad immediate basin of attraction, the boundary separating any two basins is a Lakes-of-Wada fractal: every boundary point touches all three root basins simultaneously.

At the origin z = 0 where the derivative p′(z) = 3z² vanishes, Newton's map has a pole of order 2 that slings nearby points toward infinity and wraps all three basins into six-fold symmetric chains of rosettes. Adding a complex relaxation factor a = 1 + 0.32i twists these radial necklaces into logarithmic spirals.

Newton fractal equation

The generalized Newton-Raphson iteration for a complex polynomial p(z) is:

For the canonical cubic polynomial and its derivative:

Its three attracting fixed points are the cube roots of unity:

Each pixel supplies the initial complex guess z0:

While the complex relaxation parameter a (default a = 1) controls the step size and angle:

Expanding the rational map (2z³ + 1) / (3z²) for a = 1 reduces each step to real and imaginary components:

For general cubic polynomials such as Stephen Smale's counterexample p(z) = z³ − 2z + 2, the inflection point z0 = 0 is a free critical point that lands in a superattracting period-2 cycle instead of a root:

Sub-structures

Unlike the Mandelbrot set, whose boundary folds around a critical point to reveal new structures at arbitrarily deep magnifications, a fixed Newton map is uniformly expanding on its Julia boundary and linearizes into self-similar copies of its poles and periodic cycles. The richest geometric transformations emerge when we vary the polynomial family, twist the step with complex relaxation, or introduce free critical points that capture orbits into non-convergent Smale basins.

Fractal coloring

Root convergence time: given a starting point (x0, y0), we can color each pixel by identifying which root r0, r1, or r2 the orbit converges to, the number of iterations required to enter a small ε-neighborhood of that root, and the residual distance upon arrival.

const CONVERGENCE_THRESHOLD = 1e-4;
const MAX_ITERATIONS = 1024;
const ROOTS = [
  [1, 0],
  [-0.5, Math.sqrt(3) / 2],
  [-0.5, -Math.sqrt(3) / 2]
];

function iterateUntilConvergence(x0, y0) {
  let x = x0, y = y0;
  for (let n = 1; n <= MAX_ITERATIONS; n++) {
    const r2 = x*x + y*y;
    const r4 = Math.max(1e-30, r2*r2);
    const sx = (x - (x*x - y*y) / r4) / 3;
    const sy = (y + (2*x*y) / r4) / 3;
    x -= sx;
    y -= sy;
    for (let k = 0; k < 3; k++) {
      const dx = x - ROOTS[k][0], dy = y - ROOTS[k][1];
      const d2 = dx*dx + dy*dy;
      if (d2 < CONVERGENCE_THRESHOLD) {
        return [n, d2, k, true];
      }
    }
  }
  return [MAX_ITERATIONS, Infinity, 0, false];
}

We'll set the RGBA colors of every pixel in the image based on the converged root index and the number of iterations. Shifting the palette phase by one-third of a cycle per root assigns each basin its own distinct color family.

function getColor(numIterations, rootIndex) {
  return [
    255,                // red
    numIterations * 7,  // green
    rootIndex * 85,     // blue (root basin tint)
    numIterations * 15  // alpha
  ];
}

Here's a general function for drawing Newton fractals onto canvas elements. Changing the polynomial recurrence or relaxation factor a transforms the root basins and spiral geometry.

function drawNewtonFractal(canvasEl, xRange, yRange, getColor) {
  const canvasWidth = canvasEl.width;
  const canvasHeight = canvasEl.height;
  const context = canvasEl.getContext('2d');
  const canvasImageData =
    context.createImageData(canvasWidth, canvasHeight);
  const image = canvasImageData.data;
  for (let i = 0; i < canvasHeight; i++) {
    for (let j = 0; j < canvasWidth; j++) {
      const x0 = xRange[0] +
        j*(xRange[1] - xRange[0]) / canvasWidth;
      const y0 = yRange[0] +
        i*(yRange[1] - yRange[0]) / canvasHeight;
      const [numIterations, rootDistanceSq, root, converged] =
        iterateUntilConvergence(x0, y0);
      const pixelIndex = 4*i*canvasWidth + 4*j;
      if (converged) {
        const mu = getMu(numIterations, rootDistanceSq);
        const pixels = getColor({ numIterations, mu, root });
        for (let p = 0; p < 4; p++) {
          image[pixelIndex + p] = pixels[p];
        }
      }
    }
  }
  context.putImageData(canvasImageData, 0, 0);
}

Here's the cubic Newton fractal (z³ − 1 = 0) colored by root basin and integer iteration count alone, showing discrete contour bands around each root and along the six rays meeting at the origin.

Because simple roots are superattracting fixed points where the error squares on every step (|z_n − r_k| ≈ |z_{n-1} − r_k|²), taking a double logarithm of the residual root distance produces a continuous iteration count (μ) that eliminates contour banding.

const CONVERGENCE_THRESHOLD = 1e-4;

function getMu(numIterations, rootDistanceSq) {
  const logRatio = Math.log(rootDistanceSq) /
    Math.log(CONVERGENCE_THRESHOLD);
  return numIterations - Math.log(logRatio) / Math.LN2;
}

function getColor(mu, root) {
  return [ 255, mu * 7, root * 85, mu * 15 ];
}

Here's the same Newton fractal with continuous root-convergence counts and a smooth Oklab palette. Each of the three cube roots of unity receives its own 120° Lakes of Wada phase sector in persimmon coral, celadon jade, and mineral indigo, while root-centered smoothingSteps remove contour seams.

Damped and Nova variants of Newton's method introduce a complex relaxation parameter a or an additive shift +c at each step, bridging the basins of attraction of root-finding maps with the parameter-plane islands of the Mandelbrot and Julia sets.

About

All fractals on the page are interactive and rendered when they're in view. The renderer uses WebGL by default for high performance, falling back to canvas if WebGL is not available.

Fractal image controls when hovering

• shift + scroll mouse wheel — zoom in/out
• double-click (or shift + double-click) — zoom in/out at cursor
• + / − keys — zoom in/out
• pinch — zoom in/out
• click and drag — pan the view
• two-finger swipe — pan the view on touchscreens
• one-finger swipe — scroll the page
• cmd-z — undo the last camera action
• 0 — reset view to its initial position

TODO

• Plot higher-degree polynomials (z⁴ − 1, z⁵ − 1) and Nova parameter planes
• Improve coloring code for better dynamic range
• Enable customizing colors
• Enable resetting each canvas to its original state
• Allow editing x and y ranges: text inputs or click/drag slider
• Maybe show the (x0, y0) center of each canvas
• Draw boxes when selecting parts of the canvas to zoom in on