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.
Top comments (1)
:0