Mandelbrot Set

Im not sure who will read these blogs, but I really like the material and I hope that someone will find it useful in the future. Its also a pretty good way to archive a lot of content I may want to keep in the future. So... without further adu, I guess, here is the Mandelbrot set.

The original rendering of the Mandelbrot Set by Robert Brooks and Peter Matelski. (Source) Its a little lower resolution than what we can make today.

Dispite it's name, the Mandelbrot set was first drawn in 1978 by Robert W. Brooks and Peter Matelski while exploring fractional linear transformations in the complex plane (Source). Then, shortly after, in 1980, it was rendered by Benoit Mandelbrot at the Thomas J. Watson Research Center. He was studying the parameter space of quadratic polynomials (\(ax^2 + bx + c = 0\)). Those of you who are already framiliar with the Mandelbrot set are probably able to see where this is going. In his paper titled "Fractal Aspects of the Iteration of \( z \to \lambda z (1-z) \) for complex \(\lambda\) and \(z\)" (you can read the paper here) Mandelbrot explores the iteration of a function \(f\) and how that iteration behaves with different parameters. We will essentially follow his work, but with modern tools because I am writing this on a laptop with more processing resources than IBM's best computer at the time. (They had a whole 333 megaflops in 1988!)

First of all, what does it mean to iterate a function? Well, if we have a function that takes a complex number as input and outputs another complex number $$ f : \mathbb{C} \to \mathbb{C} $$ then we can apply the function to an input \(z\in\mathbb{C}\) (read \(z\) in the set of complex numbers) and get some output \(\omega\in\mathbb{C}\). Lets look at a concrete example.


Let \(f : \mathbb{C} \to \mathbb{C}\) be a function such that $$ f(z) = 2z $$ Then we can calculate the values that different inputs get mapped to. $$ \begin{align*} f(0) &= 0 \\ f(1) &= 2 \\ f(-1) &= -2 \\ f\left(\frac{1}{2}\right) &= 02\\ f(\pi) &= 2\pi \\ f(3+4i) &= 2(3+4i) = 6 + 8i \end{align*} $$

In plain terms, a function is just a relationship that maps inputs to outputs according to some process. (That isn't the most precise but it is a good intuition for how we are thinking of functions here.) In order to iterate a funciton, we simply take apply the function to some input to obtain an output, then we apply the same function to that output to obtain another, different output. To keep track of this process, we need some notation.

If we let \(z_0 \in\mathbb{C}\) (remember, this is just saying that \(z_0\) is a complex number) then we can say that $$ \begin{align*} z_1 &= f(z_0) \\ z_2 &= f(z_1) = f(f(z_0)) \\ z_3 &= f(z_2) = f(f(f(z_0))) \\ &... \\ z_n &= f(z_{n-1}) = f(f(... f(z_0)... )) \end{align*} $$ We are just applying the same function over and over again.


Let \(f(z) = 2z\) and let \(z_0 = 1\). Then we get $$ \begin{align*} z_0 &= 1 \\ z_1 &= f(z_0) = f(1) = 2 \\ z_2 &= f(z_1) = f(2) = 4 \\ z_3 &= f(z_2) = f(4) = 8 \\ &... \\ z_n &= 2^n \end{align*} $$

Now that we have some of the basics under our belts, we can get to what the Mandelbrot set actualy is. The mandelbrot set is a set of points in the complex plane \(c\in\mathbb{C}\) such that the function \(f(z)=z^2+c\) doesn't diverge when the function is iterated starting with \(z_0=0\). When we say that something doesn't diverge, we are essentially saying that the \(z_0,z_1,z_2...\) don't "blow up" to infinity as we iterate.

For the real numbers, a simple way to check that our sequence is bounded is to find some number \(b\) such that for every \(x_i\in(x_n)\) it is true that \(x_i < b\). In other words, find a \(b\) that is bigger than every element in our sequence, then because every element is less than \(b\), they can't "blow up" to infinity. There is an analagous method for the complex numbers, where we can pick some number \(\omega\) such that for every \(z_i\in(z_n)\) it is true that \(|z_i| < |\omega|\) where \(|z|\) is the complex modulus of \(z\). (\(|a+bi|=\sqrt{a^2+b^2}\))

In order to test if a point is in the Mandelbrot set, we first need a way to generate the function we are going to iterate. We will create a function get_func(c) that will take the point of interest c and return the function we are going to iterate.

lambda z : z*z + c

Now, the most important step. The function below will take in a set of \((x,y)\) coordinates and output whether the point \(x + yi\) is in the mandelbrot set or not. This is the function that decides which pixels are colored black and which pixels are colored white in our image. It is a little complicated so lets upnpack it a bit.

  1. Convert the coordinates \((x,y)\) into a single complex number \(c=x+yi\).
  2. Create the function we are going to iterate.
  3. Initialize the variable that keeps track of the values we are iterating to produce.
  4. Iterate the function 50 times.
  5. Check to see if the outputs of the function are diverging to infinity. The way this check works is if \(|z|\) gets too large, then the \(z^2\) part of the function will start to grow faster than the \(+c\) part of the function can compensate for, causing the iteration to diverge toward infinity.

def mandelbrot(x,y):
    c = complex(x, y)       # 1
    f = lambda z : z*z + c  # 2
    v = complex(0,0)        # 3
    for _ in range(50):     # 4
        v = f(v)
        if abs(v) > 5:      # 5
            return False
    return True

Now, we have a method to determine which points are in the Mandelbrot set and which ones are not. The last thing we have to do is color the pixels in an image based on whether or not their corresponding complex number is in the Mandelbrot set or not. This chunk of code is also a bit complicated so we will break it down step by step.

  1. Define the bounds of the rectangle in the complex plane we are going to draw in our image. The Mandelbrot set is mostly centered between \(\frac{1}{2}\) and \(\frac{-3}{2}\) on the real line and \(1\) and \(-1\) on the imagenary axis.
  2. Define the resolution (in pixels) of the image we want to produce.
  3. Create the 2D array that will hold the values of the pixels. This is the array we will render at the end of the script. Each cell in the array will correspond to one pixel in the image. A 0 corresponds to a white pixel and a 1 corresponds to a black pixel.
  4. Iterate over each of the pixels \((i,j)\) in the image.
  5. Map each pixel to it's corresponding point in the complex plane. This is done by converting each pixel to a value between \(0\) and \(1\) by computing \(v=(i/WIDTH)\) (for the real axis in this example) which represets the "progess" of that pixel along the axis. Then, the point on the complex plane is computed using what is called linear interpolation. A linear interpolation is a way to smoothly transition between two values using some parameter \(0\leq v\leq 1\). The formula for interpolating between \(a\) and \(b\) is $$ va + (1-v)b $$ For example, if \(v=0\), then the result will be \(b\), if \(v=1\) then the result is \(a\). If \(v=\frac{1}{2}\) then the result is \(\frac{1}{2}a + \frac{1}{2}b\).
  6. Compute the value for that pixel and assign it in the pixel array.

import numpy as np

x_lim = (-1.5, 0.5)                                             # 1
y_lim = (-1, 1)

HEIGHT = 800                                                    # 2
WIDTH = 1200

image = np.zeros(shape=(HEIGHT, WIDTH))                         # 3

for i in range(WIDTH):                                          # 4
    for j in range(HEIGHT):
        x = (i/WIDTH)*x_lim[1] + (1-(i/WIDTH))*x_lim[0]         # 5
        y = (j/HEIGHT)*y_lim[1] + (1-(j/HEIGHT))*y_lim[0]
        if mand(x, y):                                          # 6
            image[j,i] = 1
            image[j,i] = 0

The final step is to simply plot the image so that we can see the results of our computation.

import matplotlib.pyplot as plt
plt.figure(figsize=(16, 12), dpi=80)
plt.imshow(image, cmap='Greys', vmin=0, vmax=1)


An image of the mandelbrot set in black and white.

And there it is! Our rendering of the mandelbrot set. We used a pretty simple method to determine the color of each pixel by whether it is in the set or not, but there are many other ways that you can color it.

We know that some points diverge to infinity and get arbitrarily large, but what about the other points? What do they do? When a point doesn't diverge to infinity, it actually converges to some set of values. Thats right, more than one value (sometimes)!

Convergence is a topic that should probably get its own set of posts to atequitly explain it, but we only need a bit of intuition here. For a sequence of numbers \((x_n) = \{x_1, x_2, x_3,...\} \) to converge, they should get colser and closer to some finite value. More specifically, for \((x_n)\) to converge to some value \(l\), then for every value \(\epsilon\in\mathbb{R}\) that is greater than zero, there should be an integer \(N\) also greater than zero, such that every \(x_n\) with \(n>N\), it should hold that \(|x_n - l| < \epsilon\). That is a whole mouth full of math, so to put it simply, for \((x_n)\) to converge to \(l\), the sequence should eventuially get within \(\epsilon\) of \(l\) and stay there.

Now with respect to the Mandelbrot set, some points in the set converge to a single value as described above, but there are other points that, when iterated, converge to jumping around between a few values. Maybe we can color the pixels depending on the length of the cycle of points it is converging to. Another possibility is coloring the pixels not in the set by how fast they diverge. The posibilities are actually endless