# Newton-Raphson iterative fractals.

(click on the image for a larger picture)

The Newton-Raphson iterative method for finding roots to equations is chaotic. If we use it to find the complex fourth roots of unity, i.e. solutions of the equation:

z^4 - 1 = 0

we find that for starting values close to any of the four roots (+/- 1, +/-i) the method quickly converges. However for values where y is nearly equal to +/-x, we find the method becomes extremely unpredictable, and that among these regions there exist further regions of stability.

Remember the equation is:

f(z_0)
z_1  =  z_0  -  -------
f'(z_0)
For the quartic (fourth roots), this expands to:
z = x + i y

f(z)  =  z^4 - 1

= (x + i y) ^ 4 - 1

= x^4 + 4 x^3 i y + 6 x^2 i^2 y^2 + 4 x i^3 y^3 + i^4 y^4 - 1

Re(f) = x^4 - 6 x^2 y^2 + y^4 - 1

Im(f) = 4 x^3 y - 4 x y^3

f'(z) = 4 z^3

= 4 (x + i y) ^ 3

= 4 [ x^3 + 3 x^2 i y + 3 x i^2 y^2 + i^3 y^3 ]

Re(f') = 4 [ x^3 - 3 x y^2 ]

Im(f') = 4 [ 3 x^2 y - y^3 ]
And that remains is the division of the two expressions; recall that complex division is achieved thus:
a + ib    (a + ib)(c - id)    (ac + bd) + i(bc - ad)
------  =  ---------------  = ----------------------
c + id    (c + id)(c - id)           c^2 + d^2
followed by subtraction from the original values of x and y. At each iteration, the value of z_1 is tested for proximity to each of the four roots, and a colour (or symbol) is associated with each root. The number or iterations required to get to within the required accuracy of the root can be used for shading.

Plotting graphically the root and number of iterations resulting from the method applied to each point on the Argand plane, we get a rather pretty fractal. The screen is diagonally quartered, with roots at the middle of the sides. On the diagonals is the fractal region, looking like a chain of heart-shapes. Zooming in to any part of the fractal reveals the same pattern, but transformed and distorted.

The method can be extended for any number of roots greater than two (although the computation becomes more expensive as the power increases). All roots (at any power) lie on the unit circle, but with four roots they lie at integer positions so the code is simpler.

I first wrote a program to draw this fractal at school (in BBC BASIC); at the time, a full-colour VGA image would take about 4 hours to render on an Archimedes. A compiled C version running on today's hardware takes a second or so.

I have written several versions of this program; one particularly worthy of note is this version for terminals. Start an xterm with a tiny font (e.g. -font 8x4, or alternatively select "unreadable" from the menu you get when you ctrl-right-mouse-click), resize it to fill the screen, find out the dimensions (by running resize(1)) and pass the values in using the -size option. Click here for a sample. By default it's the quartic fractal; run with -cubic for the cube roots of unity.

I've also written a little Perl script that translates the output into an XPM file, from which I've created this colour image. In fact, it only took a one-line mod to the program (to get it to output, in addition to the character for the root, one for the number of iterations too) and a change to the perl script, to get it to produce the "classic" shaded image shown at the top of this document.

I also have an MFC implementation, but since the terminal version, in conjunction with UNIX tools, generates the same image, I've not bothered to put it up here.

#### Further work

I'd like to try to find other functions (complex polynomials) that have chaotic regions under the Newton method. However a more general scheme for detecting roots is required, as in general, it will not be possible to find a function's roots analytically. Perhaps convergence can be detected at run time, and each convergence point given a colour (although this may lead to an excessive number of colours). As convergence points are found, they are looked up among all other known ones to see if they are the same.