The Mandelbrot set

Mathematical background for this notebook can be found in a separate set of beamer slides. The code below is also explained in this video.

import numpy as np
import matplotlib.pyplot as plt

Here is a function to check whether a point $c$ lies in the Mandelbrot set $M$. We put $z=0$ and then start a loop. In the loop we replace $z$ by $q_c(z)=z^2+c$, so in the $k$ th iteration we have $z=q_c^k(0)$. If at any point we find that $|z|\gt 2$, we know that $c$ is not in $M$, so we return $0$. If we still have $|z|\leq 2$ after maxiter iterations, we have high confidence (but not absolute certainty) that $c$ lies in $M$, so we return $1$.

def mandelbrot_test_a(c, maxiter=100):
    """
    Test if a point c is in the Mandelbrot set.  Return 1 if it is, 0 otherwise.
    """
    z = 0
    for k in range(maxiter):
        z = z**2 + c
        if abs(z) > 2:
            return 0
    return 1

The function above is fine for testing a single point, but we really want to check many points in parallel. If $z$ and $c$ are numpy arrays then we can calculate $z^2+c$ in parallel with no trouble. However, some entries in the array $z$ will have absolute value bigger than two and others will not. The expression abs(z) > 2 will produce a vector of true/false values rather than a single true/false value, and this prevents the if clause from working correctly. The simplest fix is just to ignore the if clause and perform the full set of maxiter iterations for all of the entries in $z$, then just check at the end for entries where $|z|\leq 2$. It will be convenient to return an array with entries 1 or 0 instead of True or False. There is a trick to convert True and False to $1$ and $0$: we just have to multiply by $1$.

def mandelbrot_test_b(c, maxiter=100):
    """
    Test if a point c is in the Mandelbrot set.  
    Return an array of zeros and ones values, where one means that the point is in the set.
    """
    z = 0
    for k in range(maxiter):
        z = z**2 + c
    return (abs(z) <= 2) * 1

The problem with the above function is that the entries in the array $z$ get very large, and will usually cause an overflow error eventually. To avoid this we need a more subtle approach. As well as the array $z$, we use another array $m$, containing zeros and ones. Initially all the entries in $m$ are equal to $1$, but at each iteration we set $m$ to (abs(z) <= 2)*1, so $m_i=0$ if we have already proved that $c_i$ is not in $M$, and $m_i=1$ if we are still unsure. At the iteration step we set $z$ to $m(z^2+c)+10(1-m)$ instead of just $z^2+c$. If we already know that $c$ is not in $M$, then $m=0$ and so the new value of $z$ is $10$. If we are still unsure whether $c$ is in $M$, then we must have $|z|\leq 2$, and the new value of $z$ will be $z^2+c$, which cannot be very big. With this arrangement there is no possibility of overflow.

def mandelbrot_test(c, maxiter=100):
    """
    Test if a point c is in the Mandelbrot set.  Return 1 if it is, 0 otherwise.
    The function maps efficiently over numpy arrays.
    """
    z = 0 * c  # array of zeros of the same shape as c
    m = z + 1  # array of ones of the same shape as c
    for k in range(maxiter):
        # Throughout this loop, m is an array of 1s and 0s and z is an array of complex numbers
        # An entry in m is 1 if the iteration starting with the corresponding c has 
        # not yet escaped from the circle of radius 2, and 0 otherwise.
        # The corresponding entry in z is the current value of the iteration if m = 1,
        # and 10 if m = 0.
        z = m * (z**2 + c) + 10 * (1 - m)
        m = (abs(z) <= 2)*1
    return m

We now set cs to be an $1000\times 1000$ matrix of evenly spaced points in $\mathbb{C}$. This is achieved by taking the broadcast sum of a row vector of real values and a column vector of imaginary values. We then apply mandelbrot_test() to get a $1000\times 1000$ array of zeros and ones, with ones in places corresponding to points in $M$. We then use plt.imshow() to display this as a picture.

n = 1000
xs = np.linspace(-2.0, 0.4, n+1).reshape(1, n+1) # row vector of real parts
ys = np.linspace(-1.2, 1.2, n+1).reshape(n+1, 1) # column vector of imaginary parts
cs = xs + 1j*ys # broadcast sum gives a matrix of complex numbers
ms = mandelbrot_test(cs, 100)
plt.imshow(ms, extent=(-2, 2, -2, 2), cmap='Purples')