Draw the Mandelbrot Set in Python

The Icon of Fractal Geometry

Even if the name is new to you, you might have seen some mesmerizing visualizations of the Mandelbrot set before. It’s a set of complex numbers, whose boundary forms a distinctive and intricate pattern when depicted on the complex plane. That pattern became arguably the most famous fractal, giving birth to fractal geometry in the late 20th century:

The discovery of the Mandelbrot set was possible thanks to technological advancement. It’s attributed to a mathematician named Benoît Mandelbrot. He worked at IBM and had access to a computer capable of what was, at the time, demanding number crunching. Today, you can explore fractals in the comfort of your home, using nothing more than Python!

Fractals are infinitely repeating patterns on different scales. While philosophers have argued for centuries about the existence of infinity, fractals do have an analogy in the real world. It’s a fairly common phenomenon occurring in nature. For example, this Romanesco cauliflower is finite but has a self-similar structure because each part of the vegetable looks like the whole, only smaller:

Self-similarity can often be defined mathematically with recursion. The Mandelbrot set isn’t perfectly self-similar as it contains slightly different copies of itself at smaller scales. Nevertheless, it can still be described by a recursive function in the complex domain.

The Boundary of Iterative Stability

Formally, the Mandelbrot set is the set of complex numbers, c, for which an infinite sequence of numbers, z₀, z₁, …, z_n, …, remains bounded. In other words, there is a limit that the magnitude of each complex number in that sequence never exceeds. The Mandelbrot sequence is given by the following recursive formula:

In plain English, to decide whether some complex number, c, belongs to the Mandelbrot set, you must feed that number to the formula above. From now on, the number c will remain constant as you iterate the sequence. The first element of the sequence, z₀, is always equal to zero. To calculate the next element, z_n+1, you’ll keep squaring the last element, z_n, and adding your initial number, c, in a feedback loop.

By observing how the resulting sequence of numbers behaves, you’ll be able to classify your complex number, c, as either a Mandelbrot set member or not. The sequence is infinite, but you must stop calculating its elements at some point. Making that choice is somewhat arbitrary and depends on your accepted level of confidence, as more elements will provide a more accurate ruling on c.

With complex numbers, you can imagine this iterative process visually in two dimensions, but you can go ahead and consider only real numbers for the sake of simplicity now. If you were to implement the above equation in Python, then it could look something like this:

>>> def z(n, c):
...     if n == 0:
...         return 0
...     else:
...         return z(n - 1, c) ** 2 + c

Your z() function returns the nth element of the sequence, which is why it expects an element’s index, n, as the first argument. The second argument, c, is a fixed number that you’re testing. This function would keep calling itself infinitely due to recursion. However, to break that chain of recursive calls, a condition checks for the base case with an immediately known solution—zero.

Try using your new function to find the first ten elements of the sequence for c = 1, and see what happens:

>>> for n in range(10):
...     print(f"z({n}) = {z(n, c=1)}")
...
z(0) = 0
z(1) = 1
z(2) = 2
z(3) = 5
z(4) = 26
z(5) = 677
z(6) = 458330
z(7) = 210066388901
z(8) = 44127887745906175987802
z(9) = 1947270476915296449559703445493848930452791205

Notice the rapid growth rate of these sequence elements. It tells you something about the membership of c = 1. Specifically, it does not belong to the Mandelbrot set, because the corresponding sequence grows without bound.

Sometimes, an iterative approach might be more efficient than a recursive one. Here’s an equivalent function that creates an infinite sequence for the specified input value, c:

>>> def sequence(c):
...     z = 0
...     while True:
...         yield z
...         z = z ** 2 + c

The sequence() function returns a generator object yielding consecutive elements of the sequence endlessly in a loop. Because it doesn’t return the corresponding element indices, you can enumerate them and stop the loop after a given number of iterations:

>>> for n, z in enumerate(sequence(c=1)):
...     print(f"z({n}) = {z}")
...     if n >= 9:
...         break
...
z(0) = 0
z(1) = 1
z(2) = 2
z(3) = 5
z(4) = 26
z(5) = 677
z(6) = 458330
z(7) = 210066388901
z(8) = 44127887745906175987802
z(9) = 1947270476915296449559703445493848930452791205

The result is the same as before, but the generator function lets you calculate the sequence elements more efficiently by using lazy evaluation. In addition to this, iteration eliminates redundant function calls for the already-computed sequence elements. As a consequence, you’re not running the risk of hitting the maximum recursion limit anymore.

Most numbers will make this sequence diverge to infinity. However, some will keep it stable by either converging the sequence to a single value or staying within a bounded range. Others will make the sequence periodically stable by cycling back and forth between the same few values. Stable and periodically stable values make up the Mandelbrot set.

For example, plugging in c = 1 makes the sequence grow without bound like you just learned, but c = -1 causes it to jump between 0 and -1 repeatedly, while c = 0 gives a sequence comprised of a single value:

It’s not obvious which numbers are stable and which aren’t, because the formula is sensitive to even the smallest change of the tested value, c. If you mark the stable numbers on the complex plane, then you’ll see the following pattern emerge:

This image has been generated by running the recursive formula up to twenty times per pixel, with each pixel representing some c value. When the magnitude of the resulting complex number was still reasonably small after all iterations, then the corresponding pixel was colored black. However, as soon as the magnitude exceeded the radius of two, then the iteration stopped and skipped the current pixel.

You may find it astonishing that a relatively simple formula that only involves addition and multiplication can produce such an elaborate structure. But that’s not all. It turns out that you can take the same formula and use it to generate infinitely many unique fractals! Do you want to see how?

The Map of Julia Sets

It’s hard to talk about the Mandelbrot set without mentioning Julia sets, which had been discovered by French mathematician Gaston Julia several decades earlier without the help of computers. Julia sets and the Mandelbrot set are closely related because you can obtain them through the same recursive formula, only with different sets of starting conditions.

While there’s only one Mandelbrot set, there are infinitely many Julia sets. So far, you always started the sequence at z₀ = 0 and systematically tested some arbitrary complex number, c, for its membership. On the other hand, to find out if a number belongs to a Julia set, you must use that number as the starting point for the sequence and pick another value for the c parameter.

Here’s a quick comparison of the formula’s terms depending on which set you’re investigating:

In the first case, c represents a potential member of the Mandelbrot set and is the only input value required because z₀ remains fixed at zero. However, each term changes its meaning when you use the formula in Julia mode. Now, c works as a parameter that determines the shape and form of an entire Julia set, while z₀ becomes your point of interest. Unlike before, the formula for a Julia set expects not one but two input values.

You can modify one of your functions defined before to make it more generic. This way, you can create infinite sequences starting at any point rather than always zero:

def sequence(c, z=0):
    while True:
        yield z
        z = z ** 2 + c

Thanks to the default argument value in the highlighted line, you can still use this function as before because z is optional. At the same time, you may change the starting point of the sequence. Perhaps you’ll get a better idea after defining wrapper functions for the Mandelbrot and Julia sets:

def mandelbrot(candidate):
    return sequence(z=0, c=candidate)

def julia(candidate, parameter):
    return sequence(z=candidate, c=parameter)

Each function returns a generator object fine-tuned to your desired starting condition. The principles for determining whether a candidate value belongs to a Julia set are similar to the Mandelbrot set that you saw earlier. In a nutshell, you must iterate the sequence and observe its behavior over time.

Benoît Mandelbrot was, in fact, studying Julia sets in his scientific research. He was particularly interested in finding those values of c that produce so-called connected Julia sets as opposed to their disconnected counterparts. The latter are known as Fatou sets and appear as dust comprised of an infinite number of pieces when visualized on the complex plane:

The image in the top-left corner represents a connected Julia set derived from c = 0.25, which belongs to the Mandelbrot set. You know that plugging a member of the Mandelbrot set into the recursive formula will produce a sequence of complex numbers that converge. The numbers converge to 0.5 in this case. However, a slight change to c can suddenly turn your Julia set into disconnected dust and make the corresponding sequence diverge to infinity.

Coincidentally, the connected Julia sets correspond to c values that generate stable sequences of the recursive formula above. So, you might say that Benoît Mandelbrot was looking for the boundary of iterative stability, or a map of all the Julia sets that would show where those sets are connected and where they are dust.

Watch how choosing different points for the c parameter on the complex plane affects the resulting Julia set:

The little moving red circle indicates the value of c. As long as it stays inside the Mandelbrot set shown on the left-hand side, the corresponding Julia set depicted to the right remains connected. Otherwise, the Julia set pops like a bubble spreading into infinitely many dusty pieces.

Did you notice how the Julia sets are changing shape? It turns out that a particular Julia set shares common visual features with the specific area of the Mandelbrot set used to seed the value of c. When you look through a magnifying glass, then both fractals will appear somewhat similar.

Okay, enough theory. Time to plot your first Mandelbrot set!

Low-Resolution Scatter Plot

A quick and dirty way to visualize the Mandelbrot set is through a scatter plot, which illustrates the relationships between paired variables. Because complex numbers are pairs of the real and imaginary components, you can untangle them into separate arrays that’ll play nicely with the scatter plot.

But first, you need to turn your Boolean mask of stability into the initial complex numbers that seeded the sequence. You can do so with the help of NumPy’s masked filtering:

16def get_members(c, num_iterations):
17    mask = is_stable(c, num_iterations)
18    return c[mask]

This function will return a one-dimensional array comprised of only those complex numbers that are stable and therefore belong to the Mandelbrot set. When you combine the functions defined so far, you’ll be able to show a scatter plot using Matplotlib. Don’t forget to add the necessary import statement at the beginning of your file:

 1import matplotlib.pyplot as plt
 2import numpy as np
 3
 4np.warnings.filterwarnings("ignore")

This brings the plotting interface to your current namespace. Now you can calculate your data and plot it:

21c = complex_matrix(-2, 0.5, -1.5, 1.5, pixel_density=21)
22members = get_members(c, num_iterations=20)
23
24plt.scatter(members.real, members.imag, color="black", marker=",", s=1)
25plt.gca().set_aspect("equal")
26plt.axis("off")
27plt.tight_layout()
28plt.show()

The call to complex_matrix() prepares a rectangular array of complex numbers in the range of -2 to 0.5 in the x-direction and between -1.5 and 1.5 in the y-direction. The subsequent call to get_members() passes through only those numbers that are members of the Mandelbrot set. Finally, plt.scatter() plots the set, and plt.show() reveals this picture:

It contains 749 points and resembles the original ASCII printout made on a dot matrix printer by Benoît Mandelbrot himself a few decades ago. You’re reliving mathematical history! Play around by adjusting the pixel density and the number of iterations to see how they affect the outcome.

High-Resolution Black-and-White Visualization

To get a more detailed visualization of the Mandelbrot set in black and white, you can keep increasing the pixel density of your scatter plot until the individual dots become hardly discernable. Alternatively, you can use Matplotlib’s plt.imshow() function with a binary colormap to plot your Boolean mask of stability.

There are only a few tweaks necessary in your existing code:

21c = complex_matrix(-2, 0.5, -1.5, 1.5, pixel_density=512)
22plt.imshow(is_stable(c, num_iterations=20), cmap="binary")
23plt.gca().set_aspect("equal")
24plt.axis("off")
25plt.tight_layout()
26plt.show()

Bump up your pixel density to a sufficiently large value, such as 512. Then, remove the call to get_members(), and replace the scatter plot with plt.imshow() to display the data as an image. If everything goes well, then you should see this picture of the Mandelbrot set:

To zoom in on a particular area of the fractal, change the bounds of the complex matrix accordingly and increase the number of iterations by a factor of ten or more. You can also experiment with different colormaps provided by Matplotlib. However, to truly unleash your inner artist, you might want to get your feet wet with Pillow, Python’s most popular imaging library.

Finding Convergent Elements of the Set

Earlier, you built a stability matrix using NumPy’s vectorization to determine which of the given complex numbers belong to the Mandelbrot set. Your function iterated over the formula a fixed number of times and returned a two-dimensional array of Boolean values. Using pure Python, you can modify this function so that it works on the individual numbers rather than a whole matrix:

>>> def is_stable(c, max_iterations):
...     z = 0
...     for _ in range(max_iterations):
...         z = z ** 2 + c
...         if abs(z) > 2:
...             return False
...     return True

It looks pretty similar to NumPy’s version from before. However, there are a few important differences. Notably, the c parameter represents a single complex number, and the function returns a scalar Boolean value. Secondly, there’s an if-condition in the loop’s body, which can terminate the iteration prematurely as soon as the magnitude of the resulting number reaches a known threshold.

You can use your new function to determine if a complex number creates a stable sequence that converges or not. That translates directly to Mandelbrot set membership. However, the result may depend on the maximum number of iterations requested:

>>> is_stable(0.26, max_iterations=20)
True
>>> is_stable(0.26, max_iterations=30)
False

For example, the number c = 0.26 is located close to the fractal’s edge, so you’ll get a wrong answer with too few iterations. Increasing the maximum number of iterations can be more accurate, revealing more details on the visualization.

There’s nothing wrong with calling a function, but wouldn’t it be nicer to leverage Python’s in and not in operators instead? Turning this code into a Python class will let you override those operators and take advantage of a much cleaner and more Pythonic syntax. Moreover, you’ll be able to retain the maximum number of iterations across many function invocations by encapsulating the state in your object.

You can take advantage of data classes to avoid having to define a custom constructor function. Here’s an equivalent class-based implementation of the same code:

# mandelbrot.py

from dataclasses import dataclass

@dataclass
class MandelbrotSet:
    max_iterations: int

    def __contains__(self, c: complex) -> bool:
        z = 0
        for _ in range(self.max_iterations):
            z = z ** 2 + c
            if abs(z) > 2:
                return False
        return True

Apart from implementing the special method .__contains__(), adding a few type hints, and moving the max_iterations parameter out from the function’s signature, the rest of the code stays the same. Assuming you saved it in a file named mandelbrot.py, you can start an interactive Python interpreter session in the same directory and import your class:

>>> from mandelbrot import MandelbrotSet
>>> mandelbrot_set = MandelbrotSet(max_iterations=30)
>>> 0.26 in mandelbrot_set
False
>>> 0.26 not in mandelbrot_set
True

Brilliant! That’s precisely the kind of information you’d need to make a black-and-white visualization of the Mandelbrot set. Later, you’ll learn about a less verbose way of drawing pixels with Pillow, but here’s a crude example for starters:

>>> from mandelbrot import MandelbrotSet
>>> mandelbrot_set = MandelbrotSet(max_iterations=20)

>>> width, height = 512, 512
>>> scale = 0.0075
>>> BLACK_AND_WHITE = "1"

>>> from PIL import Image
>>> image = Image.new(mode=BLACK_AND_WHITE, size=(width, height))