(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)
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.
