## Example A Newton Fractal

In this example we use complex math to create a fractal based on the Newton-Raphson method (or simply Newton's method). Fractals are used by scientists to investigate chaotic systems: systems whose state over time is highly dependent on initial conditions. The purpose of this example is to show the capabilities of Python's complex math and explore some ways to visualize data other than a regular plot; fractals are a perfect match.

Newton's method is an iterative procedure to find numerical solutions, or roots, to an equation of the form f(z) = 0 using an initial guess. One of the characteristics of the method is that in case of several solutions, we cannot tell in advance, based on the initial guess, what the converged solution will be (usually). If you were to map out the initial guesses based on the solution, you would find they converge to results in an image known as Newton's fractal, which is geometrically interesting.

If you'd like to read more about Newton's method, have a look at http://en.wikipedia. org/wiki/Newton%27s_method; there's a lot of additional information available on the Internet. The function we'll map is f(z) = z4 + 1. This function has four roots:

>>> solutions = [cos((2*n+l)*pi/4)+lj*sin((2*n+l)*pi/4) for n in range(4)]

To verify that these are indeed solutions to the equation: >>> [z**4 for z in solutions]

[(-l+4.4408920985006262e-0l6j), (-l+4.4408920985006262e-0l6j), (-1.0000000000000 004+6.66l338l477509412e-0l6j), (-l+8.88l7841970012523e-0l6j)]

The imaginary parts are on the order of scale of 10-16 and are due to inaccuracies of the trigonometric functions, n, and the floating-point representation; the imaginary parts are actually zero.

Newton's method takes an initial guess and calculates the next guess by applying the equation zn+1 = zn - f(zn) / f(zn), where f(z) is the derivative of f(z), or in our case z = z-(z**4+l)/(4*z**3). To check whether the new value is a "good" solution, we reapply it to the original equation, f(z), and check how close it is to zero. In reality, we check whether the absolute value is less than delta, a predefined small value. The number of iterations is an indication of how fast the solution was reached. We'll use this to select the color depth of each solution: solutions that converged fast will be brighter. Once our guess converges, we check what solution it converged to and color it accordingly. Since there are four solutions, there will be four colors (at varying color depths) in the fractal. Listing 7-1 generates said Newton's fractal in the region (0, 0)-(1, 1).

Listing 7-1. fractal.py from PIL import Image from cmath import *

# creates a z**4+l = 0 fractal using the Newton-Raphson

# root finding method delta = 0.000001 # convergence criteria res = 800 # image size iters =30 # number of iterations

# create an image to draw on, paint it black img = Image.new("RGB", (res, res), (0, 0, 0))

# these are the solutions to the equation z**4+l = 0 (Euler's formula) solutions = [cos((2*n+l)*pi/4)+lj*sin((2*n+l)*pi/4) for n in range(4)] colors = [(1, 0, 0), (0, 1, 0), (0, 0, 1), (1, 1, 0)]

for re in range(o, res): for im in range(o, res): z = (re+lj*im)/res for i in range(iters): try:

# possibly divide by zero exception continue if(abs(z**4+l) < delta): break

# color depth is a function of the number of iterations color_depth = int((iters-i)*255.0/iters)

# find to which solution this guess converged to err = [abs(z-root) for root in solutions] distances = zip(err, range(len(colors)))

# select the color associated with the solution color = [i*color_depth for i in colors[min(distances)[l]]] img.putpixel((re, im), tuple(color))

img.save('../images/fractal_z4s_%03d_%03d_%03d.png' % \ (iters, res, abs(logl0(delta))), dpi=(l50, 150))

We use the Python Imaging Library (PIL) to draw the fractal. We start by creating an RGB image of size res specifying the fractal's resolution. We then implement Newton's method with a for loop, and an if statement to check for convergence.

While the iteration is straightforward, deciding which of the four solutions a specific guess converges to and then mapping to the right color and color depth requires some clarifications.

The colors list is composed of the colors red, green, blue, and yellow, each represented by a tuple of Red-Green-Blue (RGB) values:

colors = [(l, 0, 0), (0, l, 0), (0, 0, 1), (l, l, 0)]

Variable color depth is directly responsible for the color depth (or shade) of the displayed color. For a small number of iterations, color depth is closer to 255, and for a greater number of iterations, this number is closer to 0, resulting in a brighter color for faster converging points (smaller number of iterations).

Once color depth is calculated, we find the solution closest to our converging value. Since we're using complex numbers, the value closest is the one with the minimum distance, or in complex math, the one with the smallest value of err = abs(guess-solution).

To implement this, we generate a list of values corresponding to the distances using a list comprehension. Here's an example using an arbitrary point:

(0.70710678ll88838l8+0.70710678ll88838l8j) >>> err = [abs(z-root) for root in solutions] >>> err

[3.23949326e-012, 1.41421356, 2.00000000, 1.41421356]

Next, we combine these values with the numbers 0-3, which represent the indices to the colors list, using the zip() function:

[(3.23949326e-012, 0), (1.41421356, l), (2.0000000, 2), (1.4142135, 3)]

We then find the minimum error by calling the function min(). To find the correct color, we index the color associated with min(distances)[l], which is the second element in the zipped tuple. Maybe it's easier to show interactively than explain:

>>> distances = zip(err, range(len(colors))) >>> min(distances) (3.23949326e-012, 0) >>> min(distances)[l] 0

Finally, we use a list comprehension to multiply the RGB values by the color depth. This is because the putpixel() method requires a tuple detailing the RGB values:

color = [i*color_depth for i in colors[min(distances)[l]]] img.putpixel((re, im), tuple(color))

■Tip As you experiment with parameters, you may wish to save some of the outputs. These runs can take a considerable time to complete, so it's a good idea to have different file names for the outputs as opposed to a single file name, ensuring files are not accidentally overwritten. Unlike data files, the outputs of these runs are dependent on input parameters and the code (e.g., version of the script) that generated them and are not date or time dependent. It doesn't matter whether the run was performed last week or last year; the results should be the same. The notation I've used is one that details all the parameters used to create the output within the file name: '../images/fractal_z4s_%03d_%03d_%03d.png' % (iters, res, abs(logio(delta))). Names of the output files detail the inputs that generated them. An even better approach (one that in this case will somewhat disturb the pleasing output) is annotating the images with text describing the parameters used. And lastly, if you use a version control system (see Chapter 2), the version number of the script that generated the output is a very welcomed addition either in the file name or in an annotation.

Figure 7-1 is a collage of outputs generated by the script with resolution=200 and values of iters ranging from 1 to 9 (top left is i = 1; bottom right is i = 9). We'll touch on collages in Chapter 9. Figure 7-2 is the result of a longer run with resolution=800 and iters=30. Figure 7-1. Collage of Newton's fractals with iterations from 1 (top left) to 9 (bottom right) Figure 7-2. Newton's fractal, max number of iterations equaling 30

■Tip The preceding example explores the region (0, 0)—(1, 1). If you wish to explore around the origin, that is, around (0, 0), change the line z = (re+lj*im)/res to z = ((re-res/2)+lj*(im-res/2))/res.