using Plots
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.
function mandelbrot_law(current_state, constant)
return NEXT STATE EXPRESSED IN TERMS OF THE CURRENT STATE AND THE CONSTANT
end
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.
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
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.
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.
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.
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
plot_mandelbrot(-2,2,-2,2,100)