For The Love of Fractals

I have been lurking around /r/FractalPorn so much that I finally decided to put my love for Fractals down on words. These structures are hauntingly beautiful and on this post, I want to share their beauty with you. Although not as technical as previous posts, it should be interesting nonetheless.

A fractal is a never-ending pattern. Fractals are infinitely complex patterns that are self-similar across different scales. They are created by repeating a simple process over and over in an ongoing feedback loop. Driven by recursion, fractals are images of dynamic systems – the pictures of Chaos. – fractalfoundation.org

We will go about making our own Mandelbrot Set, which is one of the easiest and most famous fractals. Here’s a beautiful example of a realtime Mandelbrot Fractal zoomer: http://jblang.github.io/XaoSjs/. Left click on an area to zoom-in and right click to zoom-out. You can see that the fractal is infinitely detailed.  Our implementation would not be as advanced; it would just display a static image.

Mandelbrot Set is a set of complex numbers c for which the function z=z² + c does not diverge. The idea is simple: we just need to perform iteration over a function. In the formula  z=z² + c, both z and c are complex numbers. We begin the iteration with z set to 0 i.e. 0+0i. c is a constant and remains the same throughout the process, whereas z‘s value is updated with every iteration.

We run the iteration for, say, a 100 times. Then, the value of z changes throughout the 100 iterations. If the value goes beyond a set limit, then we know it has diverged, whereas if it does not leave the set limit even after performing those 100 iterations, then it is said to have converged (for that 100 iterations).

Let,

z = (a +bi) and c = (d+ei)

Then,

z² = (a²-b²) + (2ab)i

z² + c =  (a²-b²) + (2ab)i + (d+e)i = (a²-b²+d) + (2ab + e)i

We now need to do the same thing that we did above for a range of complex values c. We start with a xy plane that represents the real components on the x axis and the imaginary components on the y axis. Then, we take a particular point on the plane, say (1+1i), test it for convergence, just like we did above. In order to do that, first of all it is assumed to be constant for the iteration period and then perform a iteration, say 100 times, starting with z set to 0. Then we compare the absolute value with our threshold each time and if its beyond the threshold then we say it has diverged. We then choose another point in the plane, say, (1+2i) and do the same for it . In this way we test for convergence of complex number c that is within the xy plane.

Now you should be beginning to get an idea as to how the fractal is visualized. In each of those xy values that we tested, if c converged, then we could give it one value, and if it diverged we could give it another value (read: 1 if converges, 0 if diverges or vice versa). How do we make it colorful, though?

In order to make the fractal colorful, we will store in an array the exact iteration count it took for the absolute value to grow beyond our chosen threshold. Then finally we will plot that array. Remember that any point in the plotted array is a complex number which has its real component on the x-axis, and the imaginary component on the y-axis. We’re plotting a complex number on the plane.

close all;
clear all;
clc;
xvals = linspace(-0.10, -0.90, 1000); %x axis plane size
yvals = linspace(-0.03, -0.69, 1000); %y axis plane size
xlength = length(xvals); % length of the real component
ylength = length(yvals); % length of the imaginary component
fractal = zeros(xlength, ylength); % the final image array
                        % to store the iteration count
                        % that results in color

                        
for i=1:xlength 
    for j=1:ylength % iterate through the components
        x = xvals(i); % pick a real value
        y = yvals(j); % pick an imaginary value
        complex_no = complex(x, y); % create a complex number
        fractal(i,j)=mbrot(complex_no,100);
                    % and check it for convergence/divergence
    end
end


image(fractal.') % show the fractal. the ' transposes the array
    % without it the image will display sideways
    % so its just for convenience
colormap(winter)

function itercount = mbrot(complex_no, max)
    z = complex(0,0); % initialize z = 0+0i
    for itercount=1:max % iterate till max no of iterations
        z = (z*z)+complex_no; % z squared + c
        if abs(z) > 4 % if absolute value is greater than 4
            break % exit loop else continue till you reach 
                  % the max no. of iterations
        end       % then return the iterations reached
                   % in case the z diverges, then iteration
                  % count will be less than max
                  % whereas if it does not diverage, iteration
                  % count will be equal to max.
    end
end

I have commented the code to make it as descriptive as possible. Here are some of the fractals produced from the above code. You can try it out yourself by changing the xvals, yvals and max number of iterations. For a more detailed study into Mandelbrots, you can read the book Make your own Mandelbrot by Tariq Rashid. I have shamelessly borrowed several concepts from the book. ;)

Maybe next time, we will look into how Julia sets (another fractal) relate to Mandelbrots.

PS. Don’t forget to go through this article: How Mandelbrot’s fractals changed the world