Oscillations in Newton’s fractal

complex-dynamicsfractalsnewton raphsonprogramming

I'm working on a program that draws Newton's fractal for a given polynomial. Newton's fractal is a fractal derived from Newton's root-finding method, which given some initial guess $x_0$ and function $f(x)$, the sequence

$\begin{align*}x_{n+1} = x_n – \frac{f(x_n)}{f'(x_n)}\end{align*}$

will usually converge to some root of $f$. This usually works fairly well for real numbers, but over the complex numbers it breaks, and for certain initial guesses it takes longer than average to converge or never converges. Newton's fractal is what you get when you color each initial point by which root or zero it's closest to. I wrote a simple algorithm in GLSL (pretty much C but allowing programs to run on the GPU) to draw it:

  • Take $x_0$ to be the coordinates of the current pixel.
  • Apply the Newton iteration some set number of times $N$, while storing the previous iterate.
  • If $\left | x_{n+1} – x_n \right | \leq \varepsilon$ for some small $\varepsilon$ during the iteration, break. This means the point has converged.
  • Take the final value $x_N$, and calculate its argument. In GLSL this is implemented using the function atan(y,x) which returns a value between $-\pi$ and $\pi$.
  • The argument is fed into a color ramp function. The brightness is based on the number of iterations until convergence.

The issue I'm having is with oscillations. When I try to plot the fractal for $f(z)=z^4 – 1$, one of the colors oscillates back and forth between two values. Choosing sufficiently small $\varepsilon$ does work here, but for some other polynomials like $f(z)=z^8 – 15z^4 + 16$ it does not. I'm planning on making my program draw the fractal for a user-inputted polynomial, so is there some way to fix (or at least reduce) this problem for all inputs?

Here is my code: https://www.shadertoy.com/view/DdcXDr

I've tried my best to comment it all out so it's clear what each part does.

Best Answer

It seems that for the sequences that converge to -1, some are terminated at the argument just below $\pi$ while others are terminated at the argument just above $-\pi$. It wouldn't be a problem if the variable t wasn't multiplied by $2\pi$ when you use it to set col, which appears to be your plan, but currently it is multiplied.

Related Question