Newton fractals with WebGL

It’s always exciting when a new 3Blue1Brown video emerges – Grant has an incredible sense for the intuitive visualisation of complex concepts (pun intended). His most recent on Newton fractals however really pushes the bar, and I couldn’t help but have a go at emulating it.

Before reading on, feel free to have a play with the demo here (it works best on a desktop browser) – drag the roots around the complex plane, or click the ‘auto’ button to have them move automatically. Scroll to zoom in and out.

WebGL demo

The first ever post on this blog, and the header image, is all about Newton fractals, so it feels fitting to revisit them many years later.

Without delving too deep back into the maths (see the post above and the 3b1b video for that!), a Newton fractal is so-called because it arises from the Newton-Raphson method for approximating the roots of functions – iterate the following equation, and as n \rightarrow \infty, x_n will (usually!) approach a root of the equation f(x) = 0:

\displaystyle{x_{n+1} = x_n - \frac{f(x_n)}{f'(x_n)}}

If f(x) is a polynomial, then it can be expressed as

\displaystyle{f(x) = \Pi_i (x - r_i)}

where the roots are r_i, i.e. places where f(r_i) = 0.

For, say, a cubic polynomial, the Newton-Raphson iteration step then looks like:

\displaystyle{x_{n+1} = x_n - \frac{(x - r_1)(x - r_2)(x - r_3)}{(x - r_1)(x - r_2) + (x - r_2)(x - r_3) + (x - r_1)(x - r_3)}}

Finally, to generate the fractal, the steps are:

  • Pick a complex x_0 and represent it as a pixel on the screen
  • Iterate the Newton-Raphson step many times
  • Colour the pixel depending on which root it ends up closest to

For the simple polynomial f(x) = x^3, you get the pretty fractal like this:

The circles denote the root locations – the colours are intentionally different so that they are visible.

Each root appears as a coloured circle, and the boundary between converged roots forms a complex fractal pattern.

The interesting notion I learned from the 3b1b video was the idea of modifying the polynomial by dragging the roots around the plane. Seeing the roots fly around the plane by themselves helps gain whatever intuitive understanding it is possible to extract from an inherently infinitely-detailed process.

It also makes for an excellent (if power-hungry) screensaver

While the maths is borderline magical, the code is pretty simple. To achieve interactive frame rates, the computation is all done on the GPU – in this case in a fragment shader.

The core loop is executed once per pixel, and is really simple – it’s just:

// One Newton-Raphson iteration step
vec2 iteration(vec2 z, vec2 r1, vec2 r2, vec2 r3) {
  vec2 one =   z - r1;
  vec2 two =   z - r2;
  vec2 three = z - r3;

  // 'div' and 'mul' are small complex multiplication functions
  return z - div(mul(mul(one, two), three), mul(two, three) + mul(one, three) + mul(one, two));
}

// Do lots of iterations (fewer when panning/zooming)
vec2 iterate = vec2(x, y);
for (int i = 0; i < 100; i++) {
  iterate = iteration(iterate, r1, r2, r3);
}

There is some other boilerplate to convert from pixels to complex numbers and back and assign colours to everything, but the core mathematical principle is a couple of lines long.

If you’re interested in exploring further, you can get the full code sample from GitHub:

https://github.com/jasmcole/Blog/tree/master/newton-fractals

(Disclaimer: it’s pretty messy, as it was written while a cat was climbing on the keyboard trying to eat the escape key. I tried to write a unit test to cover all possible values, but my laptop complained that it didn’t have infinite memory.)

PS – if you made it this far and enjoyed this post, you might want to come and work with me.

Leave a comment