Building fractals with the Julia programming language

#julialang #fractal #plotting

Hello everybody,

let me introduce myself first:

Highschool student.
Pretty much always had a passion to math.
Started programming around a year ago.
French (sorry if my English sounds sketchy)

As I always liked maths, I wondered how I could have fun while programming, something that I would find fun, and that I could share to others even if they don't know the (technical) details.
That leaves plenty of rooms for choice, but I end up doing lots of fractals, as the visual aspect makes it especially rewarding.

I do it in Julia, high-level, fast programming language. Be aware that this isn't a tutorial to Julia per se, a little bit on its libraries, it's just me showing off stuffs I did and how.

Packages used:

Luxor.jl
GLMakie.jl
CliffordAlgebras.jl
Lindenmayer.jl

Let's start really simple, and 2D, as a warm-up, let's just do a spiral. For that, we will plot every points on this spiral and join them. And let the user choose how many circle they want (it would be infinity to make a real fractal).

using GLMakie


function spiral_point(n)
    rot = (ℯ^(im*n)) # set up the rotation of the point
    point = rot/n
    (point.re, point.im)  # this point converge to (0, 0) as n grows
end

function spiral(n) # where n is the number of spirals
    V = Tuple{Float64, Float64}[]
    for i in 1:0.01:2π*n
        push!(V, spiral_point(i))
    end
    V
end

println("How many spirals do you want? (number in N)")
const num = parse(Int, readline()) 
spiral(num) |> lines |> display

I personally don't like that it doesn't make enough detail changes as n grows, would be cooler if it was slower first, then faster, so that it shows more details. Let's do a gradient function for that.

using GLMakie
function gradient(n)
    f(x) =  √(x/n)
end

function spiral_point(n, f)
    rot = (ℯ^(im*n))
    point = rot/f(n)
    (point.re, point.im)
end 

function spiral(n, f)
    V = Tuple{Float64, Float64}[]
    for i in 1:0.01:2π*n
        push!(V, spiral_point(i, f))
    end
    V
end

println("How many spirals do you want? (number in N)")
const num = parse(Int, readline())
spiral(num, gradient(num)) |> lines |> display

Looks much better. Useful technique for customization.

Recursion is a common technic to create fractals, let me show an example which looks really cool when n is small, and really familiar when n is big :)

using Luxor

const N = 9 # change that to change the deps
const SIZE = 2700 # change that to resize the image
const NAME = "MyLovelyFractal.png" # change that to change the name of your .png file you'll receive

function points(x, y, size, l=Tuple{Int, Int}[])
    0 == size && return x, y, [l; (x, y)]
    x, y, l = points(x, y, size-1, l)
    y -= 2^(size-1)
    x -= 2^(size-1) - 1
    x, y, l = points(x, y, size-1, l)
    x += 1
    y += 2^(size-1)
    points(x, y, size-1, l)
end
p(s) = points(0, 0, s)[end]

function main()
    v = p(NUM)
    Drawing(SIZE, SIZE, NAME)
    origin(0, SIZE)
    setcolor("blue")
    for i in 1:length(v)-1
        line(Point(v[i])*5, Point(v[i+1])*5, :stroke)
    end
    finish()
    preview()
end

main()

I'm quite proud of this little one :)
As you can see, as n grows larger, it looks like the Sierpinski triangle.

Speaking of Sierpinski triangle, and as we are in the line section, what about using Lindenmayer-System to do it?

using Lindenmayer
Sierpinski = LSystem(Dict("G" => "F+G+F", "F" => "G-F-G"), "F")

drawLSystem(Sierpinski,
    forward = 6,
    startingorientation = 0,
    turn = 60,
    iterations = 6,
    filename = "Sierpinsky.png")

For Hilbert:

Hilbert = LSystem(Dict("X" => "-YF+XFX+FY-", "Y" => "+XF-YFY-FX+"), "X")
# rotation unit = -90

and for some plants:

Plant = LSystem(Dict("F" => "FF+[+F-F-F]-[-F+F+F]"), "F")
# rotation unit = 22.5

Think about the possibilities!✨✨✨

You'll see lots of ways to draw plants when looking for fractals, so let's make another one, but this time, with random number and finally starting to scatter point-wise instead of lines (it allows for more flexibility).

using GLMakie
f1(x, y) = [0 0; 0 0.16] * [x, y]
f2(x, y) = [0.85 0.04; -0.04 0.85] * [x, y] + [0, 1.6]
f3(x, y) = [0.2 -0.26; 0.23 0.22] * [x, y] + [0, 1.6]
f4(x, y) = [-0.15 0.28; 0.26 0.24] * [x, y] + [0, 0.44]

function f(x, y)
    r = rand()
    r < 0.01 && return f1(x, y)
    r < 0.86 && return f2(x, y)
    r < 0.93 && return f3(x, y)
    f4(x, y)
end

function barnsley_fern(iter)
    X = [0.0]
    Y = [0.0]
    for i in 1:iter
        x, y = f(X[end], Y[end])
        push!(X, x)
        push!(Y, y)
    end
    scatter(X, Y, color=:green, markersize=1)
end

barnsley_fern(1_000_000)

Now that we are doing point-based fractal, we can do the mandelbrot and Julia's sets!

f(z, c) = z*z + c
function is_stable(iter, z, c)
    for _ in 1:iter
        if abs(z) > 2
            return false
        end
        z = f(z, c)
    end
    true
end

function mandel(precision, X, Y, f) # passing this f will be explained just after
    Points = Tuple{Float64, Float64}[]
    for x in X
        for y in Y
            z = f(x, y)
            if is_stable(precision, z, z)
                push!(Points, (x, y))
            end
        end
    end
    scatter(Points, markersize=1)
end


function julia(precision, X, Y, c, f)
    Points = Tuple{Float64, Float64}[]
    for x in X
        for y in Y
            z = f(x, y)
            if is_stable(precision, z, c)
                push!(Points, (x, y))
            end
        end
    end
    scatter(Points, markersize=4)
end

mandel(80, -2:0.0025:2, -2:0.0025:2, Complex) |> display
sleep(10)
julia(80, -2:0.0005:2, -2:0.0005:2, -0.8im, Complex) |> display

Julia set made in Julia, lovely.

Now... why did I pass the f? Welp, it's because you can decide to make your own to see how it will behave :)
This example is in 3D, so I didn't reuse some of the functions I made, but you will get the idea:

using CliffordAlgebras, GLMakie
const S = CliffordAlgebra(1, 1, 0)
const e1 = S.e1
const e2 = S.e2
# f definition
# is_stable definition
import Base.abs
abs(n::MultiVector) = vector(n) .^ 2 |> sum |> sqrt  # type piracy /(o_o)\

function mandel(precision, X, Y, Z, f)
    Points = Tuple{Float64, Float64, Float64}[]
    for x in X
        for y in Y
            for z in Z
                num = f(x, y, z)
                if is_stable(precision, num, num)
                    push!(Points, (x, y, z))
                end
            end
        end
    end
    Points
end
const RANGE = -2:0.05:2
fun(a, b, c) = a + b*e1 + c*e2
mandel(50, RANGE, RANGE, RANGE, fun) |> (x -> scatter(x, markersize=4))

Now you have a 3D slice of the mandelbrot set if defined in Cl(1, 1, 0).
And now, the last part, seeing this whole set thanks to an animation !

Let's suppose we have the same code until the fun(a, b, c)

n = Observable(0.0)
fn = lift(d -> (a, b, c) -> (a + b*e1 + c*e2 + d*e1*e2), n)
p = lift(f -> mandel(50, RANGE, RANGE, RANGE, f), fn)
fig, axe = scatter(p)
const frames = 1:28
record(fig, "Cl1-1-0_4DMandel.gif", frames; framerate=28) do frame
    n[] = -1.5+frame/10
end

See? with the help of passing f, you can do lots of things !
And so, here is the end, you're one of the lucky one who have seen the whole (except for all the missing points ...) mandelbrot set if defined in CliffordAlgebra(1, 1, 0) !
Hopes it was interesting.

So, why Julia?
You may not notice it because it's all compiling time, and TTFP is a real thing, but Julia is really fast, once I have my algo down, I can runs way more iterations than what I can in python.
It have some really great maths capacities, it's close to math in notation, and it feels quite natural.
It's extremely expressive, it just makes things easier to get something working well.
Awesome community, I really like some persons from the Human of Julia (discord) server.