Madelbrot fractal

The Mandelbrot set is a famous fractal that can be implmented quite easily in Python. We will use a numpy array to create the image pixels, then save the image using the technique described here.

What is the Mandelbrot set?

The Mandelbrot set is defined by the following equation:

xnext = xx - yy + c1 ynext = 2xy + c2

Starting from an initial value of x and y, and two constant values c1 and c2, we can calculate two new values, xnext and ynext.

Iteration

The next key ingredient of the Mandelbrot set is iteration. After calculating the new values of x and y, we feed them back into the equation and run the calculation again, over and over. Here is a Python loop that performs the calculation. The loop starts with x and y set to zero, and prints out the new value after each loop. We stop after 10 cycles:

def calc(c1, c2):
    x = y = 0
    for i in range(10):
        xnext = x*x - y*y + c1
        ynext = 2*x*y + c2
        x = xnext
        y = ynext
        print(x, y)

Now if we try calling this function with initial values of, for example c1 = 0.2 and c2 = 0.1, the sequence of values looks like this:

0.2 0.1
0.23 0.14
0.2333 0.16440000000000002
0.22740153000000002 0.17670904
0.22048537102861931 0.18036781212166242
0.21608125118807261 0.17953692795453008
0.21445759861565283 0.17758912805375543
0.21445416320109933 0.1761706758853121
0.21495448107239606 0.1755610697551134
0.2153837397195433 0.17547527729145027

After a few loops, the values kind of stabilise - they converge. If we try slightly different numbers, such as c1 = 0.5 and c2 = 1.0, we get something very different:

0.5 0.01
0.7499 0.02
1.06195001 0.039996000000000004
1.6261381437230003 0.09494750519992001
3.135310233727196 0.31879551971385567
10.228539678324857 2.0090457108504634
101.08675928277933 41.10920753800466
8529.065957891833 8311.20313340019
3668869.0894282013 141773599.43841496
-2.0086292897328776e+16 1040297553353152.1

The numbers get very big, very quickly - the diverge.

The Mandelbrot set

If we try this experiment with lots of different values for c1 and c2, we find that some values converge and some diverge. We could make a plot, representing each pair of values by the pixel at point (c1, c2). We could make the pixel black if the pair converges, and white if it diverges. The converging points are called the Mandelbrot set.

Since the equations are so simple, you might expect the shape created by the black pixels to also be simple - a circle or something. But in fact the shape is incredibly intricate, the edge of the shape is so complex that you can zoom in on it forever and never run out of finer and finer detail.

To plot this shape, we will first tidy up our calc function:

def calc(c1, c2):
    x = y = 0
    for i in range(1000):
        x, y = x*x - y*y + c1, 2*x*y + c2
        if x*x + y*y > 4:
            return i+1
    return 0

Here we have:

We return 0 if the loop finishes (ie if the function converges, which means the point is in the set), and none zero if the function diverges.

Plotting the image

We are going to use a numpy array to hold our image as we create it. Each element of the numpy array holds the colour value of one pixel.

It turns out that the interesting part of the Mandelbrot set is in a a region that is 3 units square, with its bottom corner at (-2, -1.5). In other words, x values between -2 and +1, and y values between -1.5 and +1.5. We would like to create an image 600 pixels square that represents this region:

scaling

The two squares represent the same region, measured in different coordinate spaces We call the first set (the 3 by 3 square) our user space, and the second set (the 600 by 600 square) our pixel space.

To create the image we need to look at every pixel, and calculate whether it is in the set or not. To do this we need to convert the pixel coordinates into user space, because these are the numbers (c1 and c2) we need to plug into our loop to do the calculation.

We define our user space using the size and start position, like this:

WIDTH = 3
HEIGHT = 3
XSTART = -2
YSTART = -1.5

To define our pixel space, we simply define a scaling factor. A PIXEL_SCALING of 200 means 3 units in user space maps onto 600 pixels in pixel space:

PIXEL_SCALE = 200

To convert a pixel position (i, j) in pixel space, into user space values c1 and c2, we do the following scaling:

c1 = XSTART + i/PIXEL_SCALE
c2 = YSTART + j/PIXEL_SCALE

The great thing about this is, if we want a more detailed image, we only have to alter the PIXEL_SCALING. For example if we change it to to 1000 our image will show the same region but at 3000 pixels resolution (but it might take a while to run).

The main loop

Now we are in a position to write our main loop. First we create our numpy array, initially filled with zeros (black pixels). Here is how we create it:

image_width = int(PIXEL_SCALE*WIDTH)
image_height = int(PIXEL_SCALE*HEIGHT)

array = np.zeros((image_height,
                  image_width,
                  3),
                 dtype=np.uint8)

The array has 3 dimensions. The first is the height, or image y axis (numpy arrays conventionally store data in row, column order rather than x, y order). The image height is pixels is the height in user space multiplied by the pixel scale (giving 600). The second dimension is the width, or x axis, which is calculated similarly. Finally, for an RGB image every pixel has 3 values (red, green and blue) so our array required a third dimension (depth) of 3. Each element of the array is an 8 bit unigned int (a byte) which cans store a value in the range 0 to 255.

Now we loop over the every pixel (i, j) in the image and do the following:

  1. Convert i and j to user space values c1 and c2
  2. Call calc to check if the point is inside the Mandelbrot set
  3. If it is not in the set, make the pixel white, otherwise leave it black

Here is the code:

for i in range(image_width):
    c1 = XSTART + i/PIXEL_SCALE
    for j in range(image_height):
        c2 = YSTART + j/PIXEL_SCALE
        v = calc(c1, c2)
        if v:
            array[j, i,] = (255, 255, 255)

Finally, we convert the numpy array to an image and store it as a PNG file (as described here).

img = Image.fromarray(array)
img.save('mandelbrot.png')

Giving this image:

scaling

Here is the full code:

import numpy as np
from PIL import Image

PIXEL_SCALE = 200
WIDTH = 3
HEIGHT = 3
XSTART = -2
YSTART = -1.5

image_width = int(PIXEL_SCALE*WIDTH)
image_height = int(PIXEL_SCALE*HEIGHT)

def calc(c1, c2):
    x = y = 0
    for i in range(1000):
        x, y = x*x - y*y + c1, 2*x*y + c2
        if x*x + y*y > 4:
            return i+1
    return 0

array = np.zeros((image_height,
                  image_width,
                  3),
                 dtype=np.uint8)

for i in range(image_width):
    c1 = XSTART + i/PIXEL_SCALE
    for j in range(image_height):
        c2 = YSTART + j/PIXEL_SCALE
        v = calc(c1, c2)
        if v:
            array[j, i,] = (255, 255, 255)

img = Image.fromarray(array)
img.save('mandelbrot.png')

Adding some colour

The calc function doesn’t just return a true or false value. It actually returns:

As you might expect, points that are a long way outside the boundary of the set tend to escape very quickly. Points that are very close to the boundary take longer to escape. So we can create a bit of extra interest in our image by setting the colour of the pixels outside the image according to the value returned by calc.

Here is a simple colour function that creates a color based on the escape number:

def create_color(v):
    values = [0, 64, 128, 196]
    b = values[v % 4] 
    g = values[(v//4) % 4] 
    r = values[(v//16) % 4]
    return (r, g, b)

As it happens, this function returns black if the value is 0, so we can replace this code in our original loop:

        if v:
            array[j, i,] = (255, 255, 255)

with this:

        array[j, i,] = create_color(v)

Giving this image:

scaling

Here is the full code:

import numpy as np
from PIL import Image

PIXEL_SCALE = 200
WIDTH = 3
HEIGHT = 3
XSTART = -2
YSTART = -1.5

image_width = int(PIXEL_SCALE*WIDTH)
image_height = int(PIXEL_SCALE*HEIGHT)

def create_color(v):
    values = [0, 64, 128, 196]
    b = values[v % 4] 
    g = values[(v//4) % 4] 
    r = values[(v//16) % 4]
    return (r, g, b)
    

def calc(c1, c2):
    x = y = 0
    for i in range(1000):
        x, y = x*x - y*y + c1, 2*x*y + c2
        if x*x + y*y > 4:
            return i+1
    return 0

array = np.zeros((image_height,
                  image_width,
                  3),
                 dtype=np.uint8)

for i in range(image_width):
    c1 = XSTART + i/PIXEL_SCALE
    for j in range(image_height):
        c2 = YSTART + j/PIXEL_SCALE
        v = calc(c1, c2)
        array[j, i,] = create_color(v)

img = Image.fromarray(array)
img.save('mandelbrot-colour.png')

Zooming in

You can zoom in on a part of the image by altering the parameters. Notice that, since WIDTH and HEIGHT are smaller, we needed to increase PIXEL_SCALE to keep the same image size (0.3*2000 is 600 pixels):

PIXEL_SCALE = 2000
WIDTH = .3
HEIGHT = .3
XSTART = -1.5
YSTART = -.15

This creates an image of the small bulb to the left of the main cardoid shape. You can see that this bulb itself has other smaller bulbs around it, and they too have their own even smaller bulbs. The detail goes on forever.

zoom

Things to try

Here are a few things to try:

For a more ambitious project, try using Tkinter, Pygame or similar to create an interactive fractal explorer where you can select parts of the fractal with the mouse and zoom in deeper.