The Mandelbrot Set

Just as Mersenne primes are one classic illustration of the use of computers in maths, so too is the Mandelbrot Set. It was discovered in 1979 by Benoit Mandelbrot, then working for IBM, and, like the logistic map it arises from the surprising behaviour of a quadratic equation. This time the equation is of a complex variable, and is

  zn+1=zn2+z0

The set is simply all points for which z does not increase towards infinity, but rather stays bounded close to the origin. It can be shown that if z ever exceeds a distance of two from the origin, it will reach infinity.

Those not liking complex numbers can instead consider the series as two coupled real series:

  xn+1=xn2-yn2+x0
  yn+1=2xnyn+y0

The sort of images it produces are well-known:

Mandelbrot Set

That took just over 30s to generate with the following code:

#!/usr/bin/python3
import numpy as np
import matplotlib.pyplot as plt

def mandel(z0,maxiter):
    i=0
    z=z0
    while (i<maxiter)and(z.real**2+z.imag**2<4):
        z=z*z+z0
        i=i+1
    return i

xmin=-2
xmax=1
ymin=-1.5
ymax=1.5
maxiter=255
size=200

M=np.ones((size,size))

img=plt.imshow(M,cmap="magma",interpolation='none',origin='lower',extent=[xmin,xmax,ymin,ymax])
img.set_clim(vmin=0, vmax=1)
plt.ion()
plt.title("(%f%+fi) to (%f%+fi)"%(xmin,ymin,xmax,ymax))
plt.show()
plt.pause(0.0001)

for i in range (0,size):
    x=xmin+i*(xmax-xmin)/size
    for j in range (0,size):
        y=ymin+j*(ymax-ymin)/size
        z=complex(x,y)
        M[j][i]=1-mandel(z,maxiter)/maxiter

    img.set_data(M)
    plt.pause(0.0001)

plt.ioff()
plt.show()

Half a minute is not fast for such a low resolution image, and it turns out the greatest amount of time is not spent generating the Mandelbrot set, but rather updating the image after every column. If this were changed to after every ten columns, then it takes about eight seconds.

    if (i%10==0):
        img.set_data(M)
        plt.pause(0.0001)

img.set_data(M)
plt.ioff()
plt.show()

That reducing the plotting frequency like this reduced the run-time by more than a factor of two clearly shows that over half the time was being spent in the plotting.

One can fiddle with the parameters at the start of the code to zoom in on regions of the Mandelbrot set.

xmin=-0.753
xmax=-0.723
ymin=0.165
ymax=0.195
maxiter=500
size=400

Detail of Mandelbrot Set

With the reduction in plotting frequency, that image also took about half a minute to produce.

The area of the Mandelbrot set is not known precisely. Theory can set upper and lower bounds, but they are quite a way apart. Can you modify the code which displays the whole set to produce an estimation of the area of the inside (the black part)? Is this likely to be an over-estimate or an under-estimate? How might it be made more accurate?

The above code is not very fast, so there is a discussion of a faster python Mandelbrot set, and a yet faster python Mandelbrot set. But speed is not everything, theory is also fun, and so there is a page on more Mandelbrot theory.