In [ ]:
using Plots

Law for the Mandelbrot Set

Our new system will be a bit more abstract: it will be just a single number, the result of evaluating a particular expression multiple times.

new_number = current_number^2 + constant

We always start with initial number equal to 0 for this system.

The constant is a parameter that will heavily influence how the numbers evolve. We will try to study this dependence.

In [ ]:
function mandelbrot_law(current_state, constant)
    return NEXT STATE EXPRESSED IN TERMS OF THE CURRENT STATE AND THE CONSTANT
end
In [ ]:
state_history = zeros(ComplexF64, 100)
constant = SET IT TO SOME VALUE BETWEEN -2 AND 2
for i in 2:100
    CALCULATE THE i-th STATE
    IF IT EXPLODES, STOP
end

plot(abs.(state_history))
xlabel!("repetition")
ylabel!("state")

As we mentioned previously, we can turn the steps we do with a computer into a function to be able to run it easier (and more often). We can turn this into a function as shown below.

In [ ]:
function mandelbrot_history(constant, length) 
    state_history = zeros(ComplexF64, length)
    for i in 2:length
        new_state = mandelbrot_law(state_history[i-1], constant)
        if abs(new_state) > 20
            state_history[i] = NaN
        else
            state_history[i] = new_state
        end
    end
    return state_history
end

Now it will be much easier to play with plots for different constants

In [ ]:
plot(abs.(mandelbrot_history(0.1,100)))
xlabel!("repetition")
ylabel!("state")

With the constant equal to zero, we do not get to move away from 0, the initial position.

With a constant that is big enough, we quickly get to a point where the evolution "explodes".

In particular, at -1 the oscillations do not decay.

At still more negative numbers, the oscillations become uneven.

Things become even more interesting with imaginary numbers.

Does it explode or not

A particularly interesting distinction appears between two classes of evolutions: for some values of the constant, the numbers explode, while for others they stay finite. It would be instructive to try to calculate after how many repetitions the number starts exploding.

In [ ]:
function mandelbrot_when_explodes(constant, length) 
    current_state = 0
    for i in 2:length
        new_state = mandelbrot_law(current_state, constant)
        if abs(new_state) > 20
            return i
        else
            current_state = new_state
        end
    end
    return length
end

We can now see that things become even more interesting when we include imaginary numbers!

However, complex numbers can not be represented on a single axis - they have both real and imaginary parts. We need a new way to look at this besides just a one-dimensional plot. We can use something called a "heatmap", a new two-dimensional type of plot. You can see these kinds of plots on the weather channel where the ground coordinates (latitude and longtitute) are the x and y axes, while the different colors represent different temperatures (or amounts of rain or something else).

Here is a neat example of how to make a two-dimensional grid of numbers (a two-dimensional list) that can then be plotted.

In [ ]:
function get_grid_of_coordinates(rmin,rmax,imin,imax)
    step = min((rmax-rmin)/500, (imax-imin)/500)
    real_coordinate = rmin:step:rmax
    imag_coordinate = imin:step:imax
    grid_of_complex_coordinates = zeros(ComplexF64, length(imag_coordinate), length(real_coordinate))
    for i_real in 1:length(real_coordinate)
        for i_imag in 1:length(imag_coordinate)
            grid_of_complex_coordinates[i_imag,i_real] = real_coordinate[i_real] + 1im * imag_coordinate[i_imag]
        end
    end
    return real_coordinate, imag_coordinate, grid_of_complex_coordinates
end

function plot_mandelbrot(rmin,rmax,imin,imax,length)
    r, i ,grid_of_complex_coordinates = get_grid_of_coordinates(rmin,rmax,imin,imax)
    heatmap(r, i, mandelbrot_when_explodes.(grid_of_complex_coordinates, length), ratio=1)
end
In [ ]:
plot_mandelbrot(-2,2,-2,2,100)