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.
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