Thick End of the Wedg

Fractals - The Mandelbrot and Julia Sets

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

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$.

I.e. with $f(z) = z^2 + c$ we evaluate $z_1=f(z_0),\ z_2 = f(z_1),\ z_3 = f(z_2),\ ...$ and test whether each $z_n$ is still inside the radius $2$ circle.

In practice, instead of evaluating $n$ from $1$ 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) : z = z^2 + c
    end
    return maxiter
end;

Packages

The Images package will be used to display the fractals.
The Colors package enables us to flexibly specify palettes.
The IndirectArrays package provides a convenient method to map to a colormap.

In [2]:
using Images, Colors, IndirectArrays

A detour on colour maps

Let's explore some colour palettes for plotting. There are a number of predefined (and named) palettes (e.g. "grays", "reds").

In [3]:
colormap("reds", 50) # defaults to 100 distinct colours but can choose your own number (here 50)
Out[3]:

You can define your own linearly changing maps in the different available colour spaces (of which there are many).

In [4]:
linspace(RGB(1,1,1), RGB(0,0,0), 80)
Out[4]:
In [5]:
linspace(HSV(0,1,1), HSV(330,1,1))
Out[5]:

A bit more complicated (but demonstrates the flexibility) - this is the way the Winston package defines a "jet" colour map.

In [6]:
jet = RGB{Float64}[
    RGB(clamp(min(4x - 1.5, -4x + 4.5), 0.0, 1.0),
        clamp(min(4x - 0.5, -4x + 3.5), 0.0, 1.0),
        clamp(min(4x + 0.5, -4x + 2.5), 0.0, 1.0))
    for x in linspace(0.0, 1.0, 80)]
Out[6]:

You can use the sequential_palette and diverging_palette functions which take as inputs the hue (in HSV colour space) of each colour as well as the number of discrete colours.

In [7]:
sequential_palette(12, 60) # 12 is the hue, 60 is the number of colours
Out[7]:
In [8]:
diverging_palette(255, 86, 10) # 255 is the left hand side hue and 86 is the right hand side hue
Out[8]:

And you can combine different colour maps into a single map by concatenation. The diverging red-grey and blue-grey ones below work quite well for the fractals so we can name and save them.

In [9]:
RedGrey = [flipdim(sequential_palette(12, 40), 1); colormap("grays", 40)]
Out[9]:
In [10]:
BlueGrey = [flipdim(colormap("blues", 40), 1); colormap("grays", 40)]
Out[10]:

Back to the fractals

First we need to create a canvas that we'll plot. Take the $2D$ complex plane and allocate to each point (or "pixel") the output of the function mandel i.e. the number of iterations, starting from that point in the plane, it takes to reach the "escape radius".

Comprehension syntax is a neat way to do this. The time taken for this calculation through an array is dependant on the number of pixels, the maximum number of iterations, and how quickly locations in the selected region "escape".

In [11]:
canvas = UInt8[mandel(complex(re, im)) for im in 1.25:-1e-3:-1.25, re in -2.25:1e-3:1.25];

Starting with a greyscale image:

In [12]:
IndirectArray(canvas, colormap("grays", 80))
Out[12]:

The greyscale image looks nice and clean. But let's add colour which can show some of the patterns a little clearer.

In [13]:
IndirectArray(canvas, BlueGrey)
Out[13]:

To control the canvas dimensions explicitly we'll define a function that takes as inputs the start (top left) and finish (bottom right) points as well as the width. The height is calculated to keep the real axis and imaginary axis on the same scale.

In [14]:
function create_ranges(start::Complex, finish::Complex, width::Int)
    incr = abs(finish.re - start.re) / (width - 1)
    reRange = start.re:incr:finish.re
    imRange = start.im:-incr:finish.im
    return reRange, imRange
end;

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

In [15]:
reRange, imRange = create_ranges(-0.210 + 0.675im, -0.170 + 0.645im, 3500)
canvas = UInt8[mandel(complex(re, im)) for im in imRange, re in reRange]
IndirectArray(canvas, BlueGrey)
Out[15]:

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 [16]:
function julia(z, c, maxiter=80)
    for n = 1:maxiter
        abs(z) > 2 ? (return n) : z = z^2 + c
    end
    return maxiter
end;

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

$c = -0.4 + 0.6i$

In [17]:
reRange, imRange = create_ranges(-1.55 + 1.15im, 1.55 - 1.15im, 3500)
canvas = UInt8[julia(complex(re, im), complex(-0.40, 0.60)) for im in imRange, re in reRange]
IndirectArray(canvas, BlueGrey)
Out[17]:

$c = 0.285 + 0.01i$

In [18]:
reRange, imRange = create_ranges(-1.55 + 1.25im, 1.55 - 1.25im, 3500)
canvas = UInt8[ julia(complex(re, im), complex(0.285, 0.010)) for im in imRange, re in reRange ]
IndirectArray(canvas, RedGrey)
Out[18]:

Animations

The Interact package enables the creation of animations in a relativeley straightforward way. These won't be interactable on the html page but you can create for yourself in Julia.

In [19]:
using Interact

We can vary the value of $c$ with "sliders". Using polar coordinates over the complex plane we could provide a fixed radius $r$ and just vary the angle $\theta$ or we can build sliders to vary both $r$ and $\theta$.

Fixed $r$, varying $\theta$

In [20]:
r = 0.7885
@manipulate for θ = slider(0:0.01:2π, value=0., label="θ")
    canvas = UInt8[julia(complex(re, im), complex(r*cos(θ), r*sin(θ))) for im=1.25:-4.5e-3:-1.25, re=-1.65:4.5e-3:1.65]
    IndirectArray(canvas, RedGrey)
end
Out[20]:

Varying both $r$ and $\theta$

In [21]:
@manipulate for θ = slider(0:0.01:2π, value=0., label="θ"), r = slider(0:0.01:2, value=0., label="r")
    canvas = UInt8[julia(complex(re, im), complex(r*cos(θ), r*sin(θ))) for im=1.25:-4.5e-3:-1.25, re=-1.65:4.5e-3:1.65]
    IndirectArray(canvas, BlueGrey)
end
Out[21]: