DEV Community

Julia-PBN
Julia-PBN

Posted on • Edited on

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
Enter fullscreen mode Exit fullscreen mode

Image description

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
Enter fullscreen mode Exit fullscreen mode

Image description

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()
Enter fullscreen mode Exit fullscreen mode

Image description

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")
Enter fullscreen mode Exit fullscreen mode

Image description

For Hilbert:

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

Image description

and for some plants:

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

Image description

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)
Enter fullscreen mode Exit fullscreen mode

Image description

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
Enter fullscreen mode Exit fullscreen mode

Image description
Image description

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))
Enter fullscreen mode Exit fullscreen mode

Image description

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
Enter fullscreen mode Exit fullscreen mode

Image description
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.

Top comments (1)

Collapse
 
uncomfyhalomacro profile image
Soc Virnyl Estela •

:0