Fractals - Mandelbrot and Julia Set

A simple demonstration of creating Mandelbrot and Julia set images using the Julia language and the definition of the sets.

Mandelbrot Set

The Mandelbrot set is defined on the complex plane as follows:

A complex number $c$ is in the Mandelbrot set if, for all positive integers $n$, $|z_n| < 2$ where $z_0 = 0$ and $z_{n+1} = z_n^2 + c$.

In practice, instead of evaluating for $n$ to $\infty$ we'll evaluate to some maximum number of iterations (defaults below to 80). If we are still testing after the maximum number of iterations we'll assume it is not included in the set.

The colours of the pixels are then defined by how many iterations it takes that point to reach the "escape radius" of 2. I.e. some escape quickly and are on one end of the colour scale, some do not escape at all (these are assumed to part of the Mandelbrot set) and are on the other end of the colour scale.

The function below is slightly modified from the main page of the Julia language website where it's used as one of the performance benchmarks.

In [1]:
function mandel(c; maxiter=80)
    z = c                                          # 1st iteration with z_0 = 0
    for n = 1:maxiter
        abs(z) > 2 ? (return n-1) : z = z^2 + c
    end
    return maxiter
end
Out[1]:
mandel (generic function with 1 method)

Use PyPlot to plot the fractals. PyPlot is a Julia interface to the Matplotlib plotting library from Python.

In [3]:
using PyPlot

There are many colour maps in matplotlib to choose from. This one works quite well.

In [4]:
get_cmap("RdGy")
Out[4]:
RdGy

Create a canvas by testing each point in a chosen array with the mandelbrot function. Comprehension syntax is a neat way to do this. The time taken to calculate through an array is dependant on the number of pixels (i.e. size of region and size of increments), the maximum number of iterations, and how quickly locations in the selected region either "escape" or terminate.

Then plot the image with imshow from PyPlot. Note that imshow puts the origin in the upper left (computer graphics convention) so the canvas is built with that convention i.e. imaginary numbers from max to min and real numbers from min to max.

In [5]:
canvas = [ uint8(mandel(complex(r,i))) for i=1.0:-1e-3:-1.0, r=-2.0:1e-3:1.0 ];
imshow(canvas, cmap="RdGy"); plt.axis("off");

To zoom in, different regions and increments can be chosen and viewed. The image below is from the "triple spiral valley" between the main bulb and the north bulb.

In [6]:
canvas = [ uint8(mandel(complex(r,i))) for i=0.675:-1e-5:0.645, r=-0.210:1e-5:-0.170 ];
imshow(canvas, cmap="RdGy"); plt.axis("off");

And zooming in quite deeply the image below is from a small region within the "seahorse valley" where there is a Mandelbrot set "island" within an embedded Julia set. The maximimu number of iterations i.e. maxiter is set to much higher to distinguish different behaviour. This increases the run time but not dramatically.

In [7]:
canvas = [ uint8(mandel(complex(r,i), maxiter=10000)) for 
          i=0.001645590:-1e-11:0.001645570, 
          r=-1.768667874:1e-11:-1.768667852 ];
imshow(canvas, cmap="RdGy"); plt.axis("off");

Julia Set

The Julia set is defined similarly to the Mandelbrot set with the difference that $c$ is now a fixed constant and the starting value of $z$ is now the pixel location. I.e.

A complex number $z_0$ is in the Julia set if, for all positive integers $n$, $|z_n| < 2$ where $z_{n+1} = z_n^2 + c$.

So each constant $c$ returns a different Julia set.

The function below follows the definition.

In [8]:
function julia(z, c; maxiter=80)
    for n = 1:maxiter
        abs(z) > 2 ? (return n-1) : z = z^2 + c
    end
    return maxiter
end
Out[8]:
julia (generic function with 1 method)

Below are a few examples for different values of $c$.

$c = -0.4 + 0.6i$

In [9]:
canvas = [ uint8(julia(complex(r,i), complex(-.40,.60))) for i=1.0:-1e-3:-1.0, r=-1.5:1e-3:1.5 ];
imshow(canvas, cmap="RdGy"); plt.axis("off");

$c = 0.285 + 0.01i$

In [10]:
canvas = [ uint8(julia(complex(r,i), complex(.285,.010))) for i=1.5:-1e-3:-1.5, r=-1.5:1e-3:1.5 ];
imshow(canvas, cmap="RdGy"); plt.axis("off");

$c = - 0.8350 - 0.2321i$

In [12]:
canvas = [ uint8(julia(complex(r,i), complex(-.8350,-.2321))) for i=1.0:-1e-3:-1.0, r=-1.5:1e-3:1.5 ];
imshow(canvas, cmap="RdGy"); plt.axis("off");

$c = - 0.8350 - 0.2321i$

In [13]:
canvas = [ uint8(julia(complex(r,i), complex(-.06,.67))) for i=1.0:-1e-3:-1.0, r=-1.5:1e-3:1.5 ];
imshow(canvas, cmap="RdGy"); plt.axis("off");

$c = -0.123 + 0.745i$. This is called Douady's rabbit fractal.

In [14]:
canvas = [ uint8(julia(complex(r,i), complex(-.123,.745))) for i=1.0:-1e-3:-1.0, r=-1.5:1e-3:1.5 ];
imshow(canvas, cmap="RdGy"); plt.axis("off");