Mandelbrot set

Definition

Let $\left( z_i \right)$ be a sequence of complex numbers defined as follow

$$ z_{i+1} = z^2_i + c $$

with $z_0 = 0$ and $c \in \mathbb C$ a fixed constant.

The Mandelbrot set is the set of all the complex numbers $c$ for which this sequence does not diverge (i.e. remains bounded in absolute value) ; in other words, the sequence tends to infinity for the numbers $c$ which do not belong to the Mandelbrot set (i.e. $\lim_{i \to +\infty}{|z_i|} = +\infty$ where $|z_i|$ is the modulus of $z_i$).

Here we will make a graphical representation of the Mandelbrot set.

Reference: Toutes les mathématiques et les bases de l'informatique, H. Stöcker, Dunod, p.696



In [ ]:

    
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import axes3d
from matplotlib import cm

EPSILON_MAX = 2.
NUM_IT_MAX = 32
Z_INIT = complex(0, 0)

def mandelbrot_version1(x, y):
    it = 0
    z = Z_INIT
    c = complex(x, y)

    # Rem: abs(z) = |z| = math.sqrt(pow(z.imag,2) + pow(z.real,2))
    while it < NUM_IT_MAX and abs(z) <= EPSILON_MAX:
        z = z**2 + c
        it += 1

    return 1 if it == NUM_IT_MAX else 0

REAL_RANGE = np.linspace(-2.0, 1.0, 800).tolist()
IMAG_RANGE = np.linspace(-1.2, 1.2, 800).tolist()

# Compute datas #############
xgrid, ygrid = np.meshgrid(REAL_RANGE, IMAG_RANGE)
data = np.array([mandelbrot_version1(x, y) for y in IMAG_RANGE for x in REAL_RANGE]).reshape(len(IMAG_RANGE), len(REAL_RANGE))

# Plot data #################
# (nice color maps: summer, magma, gist_gray, gist_yarg, gist_heat, Blues, coolwarm, copper)
fig, ax = plt.subplots(figsize=(14, 14))
ax.imshow(data, extent=[xgrid.min(), xgrid.max(), ygrid.min(), ygrid.max()], interpolation="none", cmap=cm.gray_r)

ax.set_axis_off()

# Set title and labels
plt.title("Mandelbrot set", fontsize=26)

plt.show()

A Python implementation

Note: the following Python script is also available there.

First, lets import some packages.



In [ ]:

    
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import axes3d
from matplotlib import cm

Then lets define the Mandelbrot set:



In [ ]:

    
EPSILON_MAX = 2.
NUM_IT_MAX = 32
Z_INIT = complex(0, 0)

def mandelbrot_version1(x, y):
    it = 0
    z = Z_INIT
    c = complex(x, y)

    # Rem: abs(z) = |z| = math.sqrt(pow(z.imag,2) + pow(z.real,2))
    while it < NUM_IT_MAX and abs(z) <= EPSILON_MAX:
        z = z**2 + c
        it += 1

    return 1 if it == NUM_IT_MAX else 0

def mandelbrot_version2(x, y):
    it = 0
    z = Z_INIT
    c = complex(x, y)

    # Rem: abs(z) = |z| = math.sqrt(pow(z.imag,2) + pow(z.real,2))
    while it < NUM_IT_MAX and abs(z) <= EPSILON_MAX:
        z = z**2 + c
        it += 1

    return it

mandelbrot_version1 defines the Mandelbrot set whereas mandelbrot_version2 is an alternative version that shows how fast the sequence diverges for a number $c = x + iy$.

Finally, lets display the Mandelbrot set in the complex plane (pixels are colored according to how rapidly the sequence diverges, the sooner the sequence diverges, the clearer the point).



In [ ]:

    
REAL_RANGE = np.linspace(-2.0, 1.0, 800).tolist()
IMAG_RANGE = np.linspace(-1.2, 1.2, 800).tolist()

# Compute datas #############
xgrid, ygrid = np.meshgrid(REAL_RANGE, IMAG_RANGE)
data = np.array([mandelbrot_version2(x, y) for y in IMAG_RANGE for x in REAL_RANGE]).reshape(len(IMAG_RANGE), len(REAL_RANGE))

# Plot data #################
fig, ax = plt.subplots(figsize=(14, 14))
ax.imshow(data, extent=[xgrid.min(), xgrid.max(), ygrid.min(), ygrid.max()], interpolation="bicubic", cmap=cm.Blues)

# Set title and labels
plt.title("Mandelbrot set", fontsize=26)
ax.set_xlabel("Re(c)", fontsize=18)
ax.set_ylabel("Im(c)", fontsize=18)

plt.show()

It can also be represented in 3 dimensions to illustrate the iterative process.



In [ ]:

    
REAL_RANGE = np.arange(-2.0, 1.0, 0.05).tolist()
IMAG_RANGE = np.arange(-1.2, 1.2, 0.05).tolist()

# Compute datas #############
xgrid, ygrid = np.meshgrid(REAL_RANGE, IMAG_RANGE)
data = np.array([mandelbrot_version2(x, y) for y in IMAG_RANGE for x in REAL_RANGE]).reshape(len(IMAG_RANGE), len(REAL_RANGE))

# Plot data #################
fig = plt.figure(figsize=(14, 10))
ax = axes3d.Axes3D(fig)
ax.plot_surface(xgrid, ygrid, data, cmap=cm.jet, rstride=1, cstride=1, color='b', shade=True)

# Set title and labels
plt.title("Mandelbrot set", fontsize=26)
ax.set_xlabel("Re(c)", fontsize=18)
ax.set_ylabel("Im(c)", fontsize=18)
ax.set_zlabel("Iterations", fontsize=18)

plt.show()