What will we cover in this tutorial?
In this tutorial you will learn what the Julia set is and understand how it is calculated. Also, how it translates into colorful images. In the process, we will learn how to utilize vectorization with NumPy arrays to achieve it.
Step 1: Understand the Julia set
Juila set are closely connect to the Mandelbrot set. If you are new to the Mandelbrot set, we recommend you read this tutorial before you proceed, as it will make it easier to understand.
Julia sets can be calculated for a function f. If we consider the function f_c(z) = z^2 + c, for a complex number c, then this function is used in the Mandelbrot set.
Recall the Mandelbrot set is calculated by identifying for a point c whether the function f_c(z) = z^2 + c , for which the sequence f_c(0), f_c(f_c(0)), f_c(f_c(f_c(0))), …., does not diverge.
Said differently, for each point c on the complex plane, if the sequence does not diverge, then that point is in the Mandelbrot set.
The Julia set has c fixed and and calculates the same sequence for z in the complex plane. That is, for each point z in the complex plane if the sequence f_c(0), f_c(f_c(0)), f_c(f_c(f_c(0))), …., does not diverge it is part of the Julia set.
Step 2: Pseudo code for Julia set of non-vectorization computation
The best way to understand is often to see the non-vectorization method to compute the Julia set.
As we consider the function f_c(z) = z^2 + c for our Julia set, we need to choose a complex number for c. Note, that complex number c can be set differently to get another Julia set.
Then each we can iterate over each point z in the complex plane.
c = -0.8 + i*0.34 for x in [-1, 1] do: for y in [-1, 1] do: z = x + i*y N = 0 while absolute(z) < 2 and N < MAX_ITERATIONS: z = z^2 + c set color for x,y to N
This provides beautiful color images of the Julia set.
Step 3: The vectorization computation using NumPy arrays
How does that translate into code using NumPy?
import numpy as np import matplotlib.pyplot as plt def julia_set(c=-0.4 + 0.6j, height=800, width=1000, x=0, y=0, zoom=1, max_iterations=100): # To make navigation easier we calculate these values x_width = 1.5 y_height = 1.5*height/width x_from = x - x_width/zoom x_to = x + x_width/zoom y_from = y - y_height/zoom y_to = y + y_height/zoom # Here the actual algorithm starts x = np.linspace(x_from, x_to, width).reshape((1, width)) y = np.linspace(y_from, y_to, height).reshape((height, 1)) z = x + 1j * y # Initialize z to all zero c = np.full(z.shape, c) # To keep track in which iteration the point diverged div_time = np.zeros(z.shape, dtype=int) # To keep track on which points did not converge so far m = np.full(c.shape, True, dtype=bool) for i in range(max_iterations): z[m] = z[m]**2 + c[m] m[np.abs(z) > 2] = False div_time[m] = i return div_time plt.imshow(julia_set(), cmap='magma') # plt.imshow(julia_set(x=0.125, y=0.125, zoom=10), cmap='magma') # plt.imshow(julia_set(c=-0.8j), cmap='magma') # plt.imshow(julia_set(c=-0.8+0.156j, max_iterations=512), cmap='magma') # plt.imshow(julia_set(c=-0.7269 + 0.1889j, max_iterations=256), cmap='magma') plt.show()