The Julia Set

Fractals – a subset of Euclidean space with a fractal dimension that strictly exceeds its topological dimension – are one of my fondest topics in mathematics. I have always been fascinated by fractals, and the first piece of code I ever wrote was to create fractals.

The Julia set is one of the most beautiful fractals to have been discovered. Named after the esteemed French mathematician Gaston Julia, the Julia set expresses a complex and rich set of dynamics, giving rise to a wide range of breathtaking visuals.

In this article, we will explore how to use Python to create the Julia set.

Complex Numbers

The Julia set lies on the complex plane, which is a space populated by complex numbers. Each point you see in the image above corresponds to a specific complex number on the complex plane with value: z = x + yi,
where i = √-1 is the imaginary unit.

Complex numbers have very interesting dynamics which is what gives the Julia set its complex dynamics, but for now let us take a look at real numbers, the numbers that we encounter in everyday life.

If you take the square of the real number 10, you get 100. Taking the square of 100 leads to 10,000. If we keep taking the square of the result, we get 100,000,000 in the next iteration, an astronomically large number in the next, and an even larger number in the next. Doing this enough times, we eventually tend to infinity, and we can say that the operation described is not bounded.

On the other hand, if you take the square of the real number 0.1, you get 0.01. Taking the square of 0.01 leads to 0.0001, and as with above if we keep taking the square of the result, we get 0.00000001 in the next number, a microscopically small number in the next and an even smaller number in the next. Doing this enough times we eventually tend to zero, and we can say that the operation described is bounded.

It turns out that complex numbers also behave similarly. If we take the square of the complex number 1 + i, we get 2i. If we square that we get -4. Repeating this we end up with 16, 256, 65536, 4294967296 and so on until we tend to infinity. This is similar to what happens if we iteratively square the result of the square of 10.

Likewise, if we take the square of 0.5 + 0.5i, we get 0.5i, -0.25, 0.0625, 0.00390625 and so on until we tend to 0. This is similar to what happens if we iteratively square the result of the square of 0.1.

Iterative Operations on the Complex Plane

Now instead of doing this iterative squaring operation on a single complex number, we can do this for an grid of complex numbers z lying on the complex plane. In addition to the squaring operation, we also add an arbitrary complex number c after the squaring operation to get the operation z := z + c.

After multiple iterations of z := z + c, most complex numbers in z will blow up to infinity, while some will remain bounded. During the iterative process we track how many iterations it takes for a complex number in z to blow up to infinity. It is this number of iterations measurement for each complex number in the grid that contains the image of the fractal in the image above.

Creating the Julia Set with Python

The Python code to create the Julia set is given below. The function julia_set takes the arguments c an arbitrary complex number, num_iter the number of iterations to perform, N the number of grid points on each axis in the grid, and X0 the limits of the grid.

import numpy as np
import matplotlib.pyplot as plt

def julia_set(c = -0.835 - 0.2321 * 1j, num_iter = 50, 
              N = 1000, X0 = np.array([-2, 2, -2, 2])):
    # Limits of the complex grid.
    x0 = X0[0] 
    x1 = X0[1]
    y0 = X0[2]
    y1 = X0[3]
    # Set up the complex grid. Each element in the grid
    # is a complex number x + yi.
    x, y = np.meshgrid(np.linspace(x0, x1, N), 
                       np.linspace(y0, y1, N) * 1j)
    z = x + y    
    # F keeps track of which grid points are bounded
    # even after many iterations of z := z**2 + c.
    
    F = np.zeros([N, N])    
    # Iterate through the operation z := z**2 + c.
    for j in range(num_iter):
        z = z ** 2 + c
        index = np.abs(z) < np.inf
        F[index] = F[index] + 1    
    
    return np.linspace(x0, x1, N), np.linspace(y0, y1, N), F

During each step in the for loop, after z = z ** 2 + c we check which points in z are still smaller than np.inf. The number of iterations taken by each point to blow up to infinity is recorded in F. By plotting F as a 2 dimensional image, the Julia set finally reveals itself. For example, the code below results in the image shown at the top of this post!

The light areas in the image correspond to points in the complex plane z that blow up to infinity very quickly after only several iterations, while the dark areas correspond to points that are bounded even after many iterations!

x, y, F = julia_set(c = 0.285 + 0.01 * 1j, num_iter = 200, 
                    N = 1000, X0 = np.array([-1.5, 1.5, -1.5, 1.5]))
plt.figure(figsize = (10, 10))
plt.pcolormesh(x, y, F, cmap = "gist_heat")
plt.axis('equal')
plt.axis('off')
plt.show()

As mentioned earlier, the Julia set has a rich and complex set of dynamics. This can be explored by changing the value of c! For example if we use the value c = -0.835–0.2321i, we get the visualization below, which is completely different from the previous one.

x, y, F = julia_set(c = -0.835 - 0.231 * 1j, num_iter = 200, 
                    N = 1000, X0 = np.array([-1.5, 1.5, -1.5, 1.5]))
plt.figure(figsize = (10, 10))
plt.pcolormesh(x, y, F, cmap = "gist_heat")
plt.axis('equal')
plt.axis('off')
plt.show()

In fact, we can take things one step further, and define c to be the entire complex grid, by using c = x + y instead using the value of the input argument in the code above. In this case, the resulting fractal is so famous that it has its own name — the Mandelbrot set.