Mandelbrot Set#
The Mandelbrot set is a mathematical object that is infamous for its aesthetically pleasing representations that one can find on the Internet. For example, if one opens its Wikipedia page one is greeted with the following image:
The Mandelbrot set is a very common example of a fractal: a mathematical object that is self-similar. What this means is that if you zoom in on the border of the Mandelbrot set, you will eventually find patterns that resemble the original, zoomed out object.
We will now guide you through the code you need to write in order to create your very first representation of the Mandelbrot set.
The Complex Plane#
First things first, we need to familiarise ourselves a little bit with the context in which our calculations will be taking place. The Mandelbrot set is drawn on the complex plane, which is the space where complex numbers live.
Complex numbers are mathematical entities much like real numbers or integers, except sometimes we have a harder time understanding what a complex number is because it is harder to “see it”, as opposed to, for example, the integer 3: you can pick up three pencils and realise that you are looking at a representation of the concept of three.
In mathematics, the imaginary unit, which is usually represented as \(i\), is a number that is defined to be the square root of \(-1\).
In APL, we can compute the square root of a number n with n*.5, so let us see what the square root of \(-1\) is:
¯1*.5
For a while now, Dyalog has added support for complex numbers, so we can work directly with them.
The imaginary unit is 0J1 in Dyalog, as we have seen above.
Now that we know what the imaginary unit is, we can define the complex numbers: a complex number is any number of the form \(a + bi\), where \(a\) and \(b\) are real numbers, the numbers we are used to dealing with.
In Dyalog, we write \(a + bi\) as aJb:
2J3
If, instead, you have the values of a and b in two variables:
a ← 2
b ← 3
Then you cannot write aJb directly because that looks like a variable name:
aJb
VALUE ERROR: Undefined name: aJb
aJb
∧
Instead, you either have to multiply b explicitly with 0J1 and add that to a:
a+0J1×b
Or you can use the circle primitive ○, for which some of the left arguments are useful when dealing with complex numbers.
For example, the left argument ¯11 can be used to multiply a number by the imaginary unit:
¯11○b
So, our complex number could be written as
a+¯11○b
The circle primitive can also be used to break complex numbers apart.
The part that comes before the J is called the real part of the complex number, and the part that comes after the J is called the imaginary part of the complex number.
The left arguments 9 and 11 can be used to extract the real and imaginary parts of a complex number, respectively:
⎕← c ← 2J3
9○c
11○c
Building the Set#
Now that we understand the fundamental pieces that we will be juggling in order to build the Mandelbrot set, we can now look at the basic algorithm that is used to build it. If the Mandelbrot set is a set in the complex plane, then it means the Mandelbrot set contains complex numbers in it. In order to build the set, we just need to know what is the rule that we used to decide if a number is in the set or not.
Iterative function#
That rule is really simple and it will build on this little function:
F ← {⍺+⍵*2}
F is a dyadic function that squares the right argument and then adds the left argument:
1 F 3
But F can also receive complex numbers as arguments, because complex numbers can also be added, subtracted, multiplied and divided.
(Granted, the way those operations work may seem strange at first, but the operations looking strange to you doesn’t mean they aren’t valid! If you want, take a little detour to Wikipedia’s article on complex numbers and, in particular, to the section on the operations.)
0J1 F 3
2J1 F 3
0J1 F 2J¯1
Before we tell you how we decide if a complex number is inside the Mandelbrot set, let us look at another operation that is useful.
Magnitude#
For real numbers, we are used to using | to determine the absolute value of a number:
|3 ¯3
The absolute value of a real number can also be seen as the distance between the number and 0. If we look at it like that, it becomes easy to define the “absolute value” of a complex number, which is actually its magnitude.
Let us think about the complex number 3J6, which in mathematical notation is \(3 + 6i\).
This number is 3 units to the right of 0, and then the \(+6i\) means we also go up 6 units:
3J6
/ |
/ |
/ |
/ |
/ |
/ |
/ |
0--------3
Now that we have this incredible diagram we can compute the distance between 3J6 and 0, we just have to apply the Pythagorean theorem: the distance between 3J6 and 0 is the hypotenuse of the right triangle drawn, which means we need to square the lengths of the sides, add them, and then compute the square root of that:
.5*⍨+/3 6*2
This is exactly what the magnitude | does for us:
|3J6
To Be or Not To Be#
We know what the magnitude of a complex number is, which is the final thing we needed in order to determine if a number is or is not inside the Mandelbrot set.
To determine if a number c is inside the Mandelbrot set, consider the sequence |F⍣0⍨c, |F⍣1⍨c, |F⍣2⍨c, |F⍣3⍨c, …
if the sequence of magnitudes tends to infinity,
cdoes not belong to the Mandelbrot set;if the sequence of magnitudes “stays small”, then
cbelongs to the Mandelbrot set.
Notice that the expression |F⍣n⍨c is just a compact way of writing |c F c F c F ... F c:
|c
|F⍣0⍨c
|c F c
|F⍣1⍨c
|c F c F c
|F⍣2⍨c
You get the idea, right?
Notice that, for this particular example, the sequence looks like it is growing very fast, and thus we could guess that 2J3 is not in the Mandelbrot set.
The problem is that, for numbers closer to 0, it is very, very, hard to predict if c is going to “explode” like 2J3 did or if it is going to stay close to 0.
What is more, just because you know if a given value c is in or out of the Mandelbrot set, it doesn’t mean you can predict if a number very close to c is in there or not:
c ← ¯.6853J.3
|F\60⍴c
The sequence above makes it look like c is inside the Mandelbrot set.
However, if we take a look at a number that is ever so slightly to the left of c, we get a sequence that suddenly starts growing a lot:
|F\60⍴c-0.0001
To be absolutely rigorous, we still do not know if c itself is inside the Mandelbrot set or not, we just looked at the initial values of the sequence.
Naturally, we can compute more terms of the sequence.
Let us check the 120th term in the sequence:
|F⍣120⍨c
We computed 120 terms now and it still looks like c is in the Mandelbrot set.
Surely this must mean that c is in the Mandelbrot set, right?
Wrong!
If we compute some more terms, we see that the sequence for c also “explodes” if we give it enough time:
|F⍣170⍨c
This poses a problem for us: for a given c, how can we know how many terms we should compute in order to check if c is in the Mandelbrot set or not?
The simple answer is: we can’t.
In practice, we choose a predetermined amount of iterations to perform, we compute those iterations and then look at the final computed element.
We then take that as an approximation: if the final element has already exploded, then c is not in the Mandelbrot set.
If the final element has not exploded, we take that c to be in the Mandelbrot set.
Sampling the Plane#
We have defined all we are going to work with, so now we can start writing the code for our first pictorial representation of the Mandelbrot set. The first step we need to take is sample the complex plane: this means that we need to build a grid of complex numbers and then we will figure out which of the points in the grid are and aren’t in the Mandelbrot set.
Building the Grid#
We start this process by defining the length of the grid, in number of points to be considered:
n ← 5
We start with a small value for n because we are still figuring out the code we will be writing, and in case we make a mistake, it is preferable that the arrays we are dealing with are small – this makes it easier to debug our program.
When we are confident in our program we can increment the value of n to create a larger grid.
After defining the size of the grid, we need to create it.
We will be creating a grid in a square that has unit length, so we can start by splitting the side of the square in n parts:
⎕← ticks ← n÷⍨¯1+⍳n+1
Think of these numbers as ticks in the horizontal axis. The plane we are working in is the complex plane, so the horizontal axis represents the real parts of the numbers in our grid and the vertical axis will represent the imaginary parts of the numbers in our grid.
Because we already have a set of ticks we can also use it for the ticks along the vertical axis. We just have to combine the horizontal and vertical ticks, to create all complex numbers that have a real part in the horizontal ticks and an imaginary part in the vertical ticks.
The vertical ticks, corresponding to the imaginary parts, are simply the horizontal ticks multiplied by the imaginary unit:
0J1×ticks
Now we just need to combine all the horizontal ticks with the vertical ticks, but we can do so by using the outer product operator:
ticks∘.+0J1×ticks
The only thing that is slightly off is that the real parts are varying vertically and the imaginary parts are varying horizontally, but we can easily fix that with the transpose ⍉:
⍉ticks∘.+0J1×ticks
Finally, if we want the grid to really mimic the complex plane, we need the upper rows of the grid to be the ones with the larger imaginary values, and not the other way around. As it stands right now, the first row of the grid has the numbers with imaginary part equal to 0 and the last row has the numbers with imaginary part equal to 1. To fix this, all we have to do is use the reverse first function:
⎕← square ← ⊖⍉ticks∘.+0J1×ticks
Adjusting the grid#
The next step lies in translating the grid and resizing it, so that it covers the area in which we know the Mandelbrot set lies.
This is relevant because we know that the Mandelbrot set lies within a square that spans from ¯2.1 and 0.9 in the horizontal axis, and ¯1.5 and 1.5 in the vertical axis.
In other words, we want the lower left corner of our grid to be ¯2.1J¯1.5 and the upper right corner to be 0.9J1.5.
Another possibility is to say that we want the length of the side of the grid to be 3 and the centre to be the point ¯0.6J0.
Our grid currently has side length 1, so we can fix the side length by multiplying by the final side length:
3×square
Finally, we have to fix the centre to be ¯0.6J0.
The current centre of the grid is 1.5J1.5, so what we can do is subtract its current centre and then add the new one:
¯0.6J0+¯1.5J¯1.5+3×square
Another possibility is recognising that subtracting a centre and adding a new one is the same as adding their difference:
¯0.6J0-1.5J1.5
We can now use their difference to shift the centre of the grid in a single operation:
⎕← grid ← ¯2.1J¯1.5+3×square
Compute the Mandelbrot Set on the Grid#
The next thing we need to do now is apply our function F to the whole grid.
Luckily, APL’s array processing capabilities allow us to do so very easily:
⎕← magnitudes ← |F⍣9⍨grid
In the matrix above we can see that we already have some really large numbers, and that is why we only apply F 9 times.
After computing applying the iterative method, we need to figure out which values are inside the Mandelbrot set and which aren’t, so we can just use a comparison:
2>magnitudes
The 1s above show points inside the Mandelbrot set and the 0s show points outside.
To make a pictorial representation of the set, we could use indexing:
' #'[1+2>magnitudes]
The representation above doesn’t look very interesting, but bear with me for a second more.
Recap#
Let us take a step back and look at the code we have written so far:
we created ticks in order to prepare for the grid:
n÷⍨¯1+⍳n+1
we then used those ticks to create a grid with complex numbers:
t∘.+0J1×t←n÷⍨¯1+⍳n+1
the grid was then readjusted, both in position and in alignment:
¯2.1J¯1.5+3×⊖⍉t∘.+0J1×t←n÷⍨¯1+⍳n+1
after that, we applied our iterative function and found out the magnitudes of the values we computed:
|{⍺+⍵*2}⍣9⍨¯2.1J¯1.5+3×⊖⍉t∘.+0J1×t←n÷⍨¯1+⍳n+1
finally, we used a comparison to turn those magnitudes into a pictorial representation of the Mandelbrot set:
' #'[1+2>|{⍺+⍵*2}⍣9⍨¯2.1J¯1.5+3×⊖⍉t∘.+0J1×t←n÷⍨¯1+⍳n+1]
If the image doesn’t look interesting, it is (in part) because our n was really small!
Now that it looks like our code is working, why don’t you make it larger?
n ← 50
⎕← M ← ' #'[1+2>|{⍺+⍵*2}⍣9⍨¯2.1J¯1.5+3×⊖⍉t∘.+0J1×t←n÷⍨¯1+⍳n+1]
Doesn’t that look much nicer? I am not saying it competes with the Wikipedia image from the beginning, but take a look at the amount of code we wrote in order to get this picture! We barely broke a sweat!
One of the things that looks odd about the representation above is its aspect ratio, because characters are generally taller than they are wide. Duplicating each character horizontally can improve that slightly:
2/M
Next Steps#
Now that you computed your very first Mandelbrot set, what should you do? There are plenty of things you could try and I will leave you with some suggestions.
Fiddle With the Numbers#
The code we wrote creates a square grid of side length 3. What happens if you change the side length? What if you make it a rectangular grid instead of a square one? What if you pick a different centre?
Create Better Pictorial Representations#
In order to be able to create more interesting representations of the Mandelbrot set, it is common to follow more or less the same workflow:
define a threshold,
t ← 2for example;define a maximum number of iterations,
iter ← 9for example;apply the iterative function a maximum of
itertimes to the grid;determine what was the iteration number for which the magnitude of the point first exceeded the threshold;
do something with the information above.
Now, instead of only caring about whether or not a point blew up during the calculation of the sequence, we also care about how quickly each point blew up.
Using the small grid from before, here is what this could look like in code:
t ← 2
iter ← 9
grid
⎕← mat ← ⊃iter-+/2<|F\9⍴⊂grid
The matrix above gives you the number of the iteration for which that specific point was such that its magnitude was above t.
This means that a very small integer value means that that point blew up very quickly, whereas a larger integer value means that that point took longer to blow up.
The points for which their value equal iter were the ones whose magnitude never surpassed the threshold:
iter=mat
This matches our first uninteresting pictorial representation of the Mandelbrot set:
' #'[1+iter=mat]
Character-Based Representations#
If you get creative with Unicode characters, you can produce the following:
iters ← iter-⊃+/2<|F\iter⍴⊂¯2.1J¯1.5+3×⊖⍉t∘.+0J1×t←n÷⍨¯1+⍳n+1
' ░▒▓█'[(4÷⍨¯1+⍳5)⍸iters÷iter]
Can you see what we are doing there? We first calculated how much it takes each point to blow up. Here is a section of those results:
⎕← section ← iters[20+⍳5;10+⍳5]
Then we need to transform the integers in the range 0 1 … iter into numbers that can be used to index into the character vector ' ░▒▓█', so we normalise the iteration counts to be decimal values between 0 and 1:
section÷iter
At the same time, we split the interval \([0, 1]\) in 4 equal segments:
(4÷⍨¯1+⍳5)
Now we can use the interval index function to map the decimal values from before into an appropriate index for ' ░▒▓█', with the caveat that only points that never blew up get mapped to █:
' ░▒▓█'[(4÷⍨¯1+⍳5)⍸section÷iter]
Colour-Based Representations#
In order to produce beautiful images like the one in the beginning of the notebook, we need to be able to map each point to a colour instead of to a character.
There are many ways to transform these integer values into colours. The simplest one would be to assign each point to a shade of grey:
⌊255×iter÷⍨iter-mat
The matrix above would work if you were using integers from 0 to 255 (like in the RGB colour scheme) but you could also use values from 0 to 1:
iter÷⍨iter-mat
You would then want to look at, for example, the documentation for Bitmap (available under Dyalog’s documentation, “Object Reference” ⇨ “(A-Z)” ⇨ “Bitmap Object”) to see how to create actual images and save them to files.
Julia Set#
There is a family of fractals that are very similar, in nature, to the Mandelbrot set, and those are called Julia sets.
In order to draw a Julia set, you do the exact same thing that we have done for the Mandelbrot set, but instead of using the function F like we have seen, you bind it to a left argument:
⎕← M ← ' #'[1+2>|0J¯0.8∘F⍣9⊢¯2.1J¯1.5+3×⊖⍉t∘.+0J1×t←n÷⍨¯1+⍳n+1]
You could explore different left arguments to bind to F or perhaps look at some of the values listed under the “Examples” section on Wikipedia.
Another thing you could do is, can you redraw this Julia set with the block characters ' ░▒▓█'?
It is not necessarily hard, but you need to tweak the respective code that works for the Mandelbrot set.
Refine the Estimation#
In our examples we applied our function F only 9 times.
Given the grid that we used, if we apply it just one more time we get a DOMAIN ERROR:
⎕← ' #'[1+2>|{⍺+⍵*2}⍣10⍨¯2.1J¯1.5+3×⊖⍉t∘.+0J1×t←n÷⍨¯1+⍳n+1]
DOMAIN ERROR
⎕←' #'[1+2>|{⍺+⍵*2}⍣10⍨¯2.1J¯1.5+3×⊖⍉t∘.+0J1×t←n÷⍨¯1+⍳n+1]
∧
The problem is that some of the numbers along the edge of the grid are growing really fast, and at some point the values become too large for APL to represent them:
⊃grid
F⍣9⍨⊃grid
F⍣10⍨⊃grid
DOMAIN ERROR
F[0] F←{⍺+⍵*2}
∧
When creating coloured representations of these fractals, you generally want the gradient of the image to be smooth.
For that to happen, you need to be able to compute more terms of the iterative sequence, so you need to come up with a way to circumvent these DOMAIN ERRORs.
A possibility is to stop applying the iterative process to points whose magnitude is already too large.
How can you do that?
Another possibility is to rewrite F slightly so that it doesn’t perform the multiplications with the large numbers.
In order to do that, you could mask out the numbers that are already too large:
F ← {m ← 2>|⍵ ⋄ (⍺×~m)+m×⍺+2*⍨⍵×m}
iter ← 20
iters ← iter-⊃+/2<|F\iter⍴⊂¯2.1J¯1.5+3×⊖⍉t∘.+0J1×t←n÷⍨¯1+⍳n+1
' ░▒▓█'[(4÷⍨¯1+⍳5)⍸iters÷iter]