Mandelbrot and Fractal
The Mandelbrot set, named after its discoverer, the French mathematician Mandelbrot, is a fractal, an infinitely ramified mathematical object that contains structure within structure within structure, as deep as we care to look. The definition of the Mandelbrot set is in terms of complex numbers as follows. Consider the equation
z' = z2 + c
where z is a complex number and c is a complex constant. For any given value of c this equation turns an input number z into an output number z′.
The definition of the Mandelbrot set involves the repeated iteration of this equation: we take an initial starting value of z and feed it into the equation to get a new value z′.
Then we take that value and feed it in again to get another value, and so forth.
The Mandelbrot set is the set of points in the complex plane that satisfies the following definition:
For a given complex value of c, start with z = 0 and iterate repeatedly. If the magnitude |z| of the resulting value is ever greater than 2, then the point in the complex plane at position c is not in the Mandelbrot set, otherwise it is in the set.
In order to use this definition one would, in principle, have to iterate infinitely many times to prove that a point is in the Mandelbrot set, since a point is in the set only if the iteration never passes |z| = 2 ever. In practice, however, one usually just performs some large number of iterations, say 100, and if |z| hasn’t exceeded 2 by that point then we call that good enough.
Write a program to make an image of the Mandelbrot set by performing the iteration for all values of c = x + iy on an N x N grid spanning the region where −2 ≤ x ≤ 2 and −2 ≤ y ≤ 2. Make a density plot in which grid points inside the Mandelbrot set are colored black and those outside are colored white.
You will probably find it useful to start off with quite a coarse grid, i.e., with a small value (perhaps N = 100) —so that your program runs quickly while you are testing it. Once you are sure it is working correctly, increase the value of N to produce a final high-quality image of the shape of the set.
The simplest way is to write it directly, we just need to create an two dimensional array to go over the given space,
and perform an iteration over z=0 via
z' = z2 + c.
If the magnitude of
from numpy import linspace, zeros, absolute
from matplotlib import pyplot as plt
a1,b1,points1 = -2,0.5,300 # Side of x
a2,b2,points2 = -1,1,300 # Side of y
N_iter = 50 # the maximum number of iterations
x_set = linspace(a1,b1,points1)
y_set = linspace(a2,b2,points2)
# Make an array to store the grids
mandset = zeros([points1,points2],float)
# Calculate the values in the array
for i in range(len(y_set)):
for j in range(len(x_set)):
#initial conditions
z = 0
c = x_set[j] + y_set[i]*1j
#we need to iterate it to see if the point satisfies the condition
for k in range(N_iter):
z = z**2+c
#print('values: ', k,z,c)
if absolute(z)>2:
mandset[i,j]=k
break
# Make the plot
plt.imshow(mandset,origin="lower",extent=[a1,b1,a2,b2],cmap="hot")
#plt.gray()
plt.show()
#plt.savefig('mandelbrot.pdf',format='pdf')
A beauty of the Mandelbrot set lies in that you can find a lot of details if you zoom in to a particular region. A high resolution is required in this case. We will also need to make the plot larger in order to see more details.
from numpy import linspace, zeros, absolute
from matplotlib import pyplot as plt
import timeit #check the total time spent
start = timeit.default_timer() #time_start
a1,b1,points1 = -2,2,1000 # Side of x
a2,b2,points2 = -2,2,1000 # Side of y
N_iter = 50 # the maximum number of iterations
x_set = linspace(a1,b1,points1)
y_set = linspace(a2,b2,points2)
# Make an array to store the grids
mandset = zeros([points1,points2],float)
# Calculate the values in the array
for i in range(len(y_set)):
for j in range(len(x_set)):
#initial conditions
z = 0
c = x_set[j] + y_set[i]*1j
#we need to iterate it to see if the point satisfies the condition
for k in range(N_iter):
z = z**2+c
#print('values: ', k,z,c)
if absolute(z)>2:
mandset[i,j]=k
break
# Make the plot
fig = plt.figure(figsize=(10, 10), dpi=150) # make the plot larger
ax = fig.add_axes([0.0, 0.0, 1.0, 1.0], frameon=False, aspect=1)
plt.imshow(mandset,origin="lower",extent=[a1,b1,a2,b2],cmap="hot")
#plt.gray()
plt.show()
#plt.savefig('mandelbrot.pdf',format='pdf')
stop = timeit.default_timer() #time_end
print('Time elapsed: ', stop - start, 'seconds')
In this code, note that we increased the number of points to 1000 * 1000 and show the plot with larger size. Obviously, it makes the plot much more attractive. To make sure the code is sufficiently fast, we also used the timeit function to track the total run time spent for the calculation.
For this simple run, it returns:
Time elapsed: 17.8881178519008 seconds
The above code is very straight forward to implement. However, it won't be very efficient if you have many points to calculate. Suppose you have N*N points to calculate with the maximum iteration of 100 times, you will need to run N*N * 100 times. This is basically a 3-for-loop structure.
As the high level language, python is not really efficient to deal with large loops. So a better way is to vectorize the code with the help of numpy. Below is an example from Joshua Adams in my class. The original code could be found in his wiki page.
import numpy as np
import matplotlib.pyplot as plt
from matplotlib import colors
import timeit
xmin, xmax, xn = -2.0, +2.0, 1000
ymin, ymax, yn = -2.0, +2.0, 1000
maxN = 200
horizon = 2
def mandelbrot_set(xmin, xmax, ymin, ymax, xn, yn, maxN, horizon=2.0):
X = np.linspace(xmin, xmax, xn, dtype=np.float32)
Y = np.linspace(ymin, ymax, yn, dtype=np.float32)
C = X + Y[:, None]*1j
N = np.zeros(C.shape, dtype=int)
Z = np.zeros(C.shape, np.complex64)
for n in range(maxN):
I = np.less(abs(Z), horizon)
N[I] = n
Z[I] = Z[I]**2 + C[I]
N[N == maxN-1] = 0
return Z, N
start = timeit.default_timer()
log_horizon = np.log(np.log(horizon))/np.log(2)
Z, N = mandelbrot_set(xmin, xmax, ymin, ymax, xn, yn, maxN, horizon)
with np.errstate(invalid='ignore'):
M = np.nan_to_num(N + 1 - np.log(np.log(abs(Z)))/np.log(2) + log_horizon)
dpi = 150
width = 10
height = 10
fig = plt.figure(figsize=(width, height), dpi=dpi)
ax = fig.add_axes([0.0, 0.0, 1.0, 1.0], frameon=False, aspect=1)
light = colors.LightSource(azdeg=90, altdeg=90)
M = light.shade(M, cmap=plt.cm.hot, vert_exag=1.5,
norm=colors.PowerNorm(.2))
plt.imshow(M, extent=[xmin, xmax, ymin, ymax])
ax.set_xticks([])
ax.set_yticks([])
plt.show()
stop = timeit.default_timer()
print('Time elapsed: ', stop - start, 'seconds')
He managed to use only 1-for-loop to get the code run. For this simple run, it returns:
Time elapsed: 2.532438270698094 seconds
Thanks to the vectorization, it only took a couple of seconds to get a 1000*1000 resolution. Well done!
from numba import jit
import numpy as np
import matplotlib.pyplot as plt
#to suppress some numba warnings
import warnings
warnings.filterwarnings("ignore")
@jit
def mandelbrot(RE ,IM, MaxIter):
c = complex(RE, IM)
z = 0.0j
for i in range(MaxIter):
z = z*z + c
if ((z.real**2 + z.imag**2) >= 4):
return i
return MaxIter
@jit
def plotTheBrot():
columns = 2000
rows = 2000
results = np.zeros([rows , columns])
for row, RE in enumerate(np.linspace(-2 , 1, rows)):
for col, IM in enumerate(np.linspace(-1, 1, columns)):
results[row , col] = mandelbrot(RE, IM, 1000)
plt.figure(dpi = 100)
plt.imshow(results.T, cmap='binary', interpolation= 'bilinear', extent=[-2,1,-1,1])
plt.xlabel('RE')
plt.ylabel('IM')
%timeit plotTheBrot()
4.21 s ± 58.5 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)