As a young person, I was deeply influenced by the book Chaos: Making a New Science by James Gleick that my father purchased. The book presents the principles behind chaos theory in an accessible way and includes stunning images of various types of fractals. I vividly recall being fascinated by the images of the Mandelbrot set and the repetitive colorful patterns. However, what truly left a lasting impression on me was the unexpected behavior of a simple mathematical expression – the logistic map.
xn+1=rxn(1−xn)
As can be seen, the map is essentially a sequence where each element depends on the value of the previous element (xn) and on the parameter r. Each new element in the sequence is the product of the previous element (xn) multiplied by r and by the value of the difference between 1 and the previous element (1−xn). The sequence contains only elements between 0 and 1. In order to ensure that the sequence stays within this range, r must be between 0 and 4.
Reproduction and starvation
The book relates to the expression as a mathematical description of an ecological system. Let’s take rabbits for example. Each element in the sequence describes the population size of rabbits in each generation. x0 is the initial population size, x1 is the population size after the first generation, x2 is the population size after the second generation, and so on. The rabbit population in the system is described by a number between 0 (indicating extinction) and 1 (indicating the largest possible population size).
On the one hand, in each generation, the number of rabbits increases at their natural reproduction rate. This is represented by the r coefficient. For example, if the number of rabbits doubles each generation (after all, rabbits…), r would be 2.
On the other hand, in each generation, there is natural rabbit mortality depending on the number of rabbits in the previous generation. There exists an inverse relationship between the number of rabbits and the natural mortality rate. The more rabbits there are, the higher the natural mortality rate due to density, competition for food, or the abundance of predators. Conversely, the fewer rabbits there are, the lower the mortality rate because there is more space and food for the few. This relationship is represented by the (1−xn) expression, which is smaller the larger xn is, and vice versa.
Let’s use some Python
It turns out that for different values of r, the logistic map sequence behaves differently. This Python program illustrates graphically the behavior of the sequence for a given r. The initial value of the sequence is not particularly important, and arbitrarily, the sequence begins with x0=0.5.
from turtle import *
from numpy import arange
r = 3.1
setworldcoordinates(0,0,20,1.1)
ht(); tracer(10000)
for y in arange(0,1.1,0.1):
pu(); setpos(0,y)
write(round(y,1))
pd(); setpos(20,y)
pu()
width(2); color("red")
x = 0.5
for n in range(100):
setpos(n,x); pd()
x = r*x*(1-x)
update()
done()For r between 0 and 1, the limit of the sequence is 0. Here is, for example, the progression of the sequence when r=0.5:
This is not surprising since the element xn is multiplied again and again by two numbers smaller than 1 (r and 1-xn).
For r values between 1 and 3, the sequence stabilizes at a constant value. For example, when r=2.3 the plot looks like this:
When r crosses the value 3, the sequence oscillates between two values. For example, r=3.1:
Attractors
Then, as r increases gradually, the sequence oscillates between 4 values, then 8 values, 16 values, and so on. These values are knows as attractors. This is with r=3.5 that has 4 attractors:
Strange attractors
However, beyond the point of r=3.57, chaos suddenly emerges! The sequence stubbornly refuses to stabilize on any values clearly, and the number of attractors gets very large with no clear order. The plots above show the progression of the different sequences between x0 to x20. Here is the plot of r=3.9 between x0 to x100:
Bifurcation diagram
Here is a python code that draws a diagram showing the converged values obtained for r. In mathematics this is called a bifurcation diagram. For each r (on the horizontal axis), the code marks vertical points at the observed attractors. For example, for r=3.1, it will mark one point at a height of 0.558 and another point at a height of 0.765, as r=3.1 has two attractors as shown before.
from turtle import *
from time import *
from numpy import arange
def draw(minr,max): # drawing the bifurcation diagram
clear()
margin=(maxr-minr)*0.02
setworldcoordinates(minr-margin,-0.05,maxr+margin,1.05) # set world coordinates with some margins
pu(); ht(); tracer(10000)
for x in arange(0,1.1,0.1): # drawing axis
setpos(minr-margin,x) #
write(round(x,1)) #
for r in arange(0,4.2,0.2): #
setpos(r,-0.02) #
write(round(r,1),align='center') #
setpos(minr,-0.04) #
write(round(minr,6),align='center') #
setpos(maxr,-0.04) #
write(round(maxr,6),align='center') #
update() #
pixel = (maxr-minr)/window_width() # calculate pixel width
for r in arange(minr,maxr,pixel): # go over all the r's in the range
x = 0.5 # x0=0.5 arbitrarily
elements = [] # elements has all the values already visited in the series
converged = False
while True:
x = r*x*(1-x) # the logistic map
if converged: # if the sequence has converged...
setpos(r,x); dot(1) # ...set a point at the relevant attractor
if any(abs(x-xn)<0.000001 for xn in elements): # check if this is a repeating element
if converged: # if the sequence has already converged...
break # ...exit
else: # if the sequence hasn't converged yet...
converged = True # ...now it is converged
elements = [x] # ...restart tracking the elements
else: # if this is not a repeating element...
elements.append(x) # ...just add it to the list
first_click = 0 # zoom in logic
def click(x,y): #
global first_click #
if x<0 or x>4: #
retrun #
pu(); setpos(x,0) #
pd(); setpos(x,1); pu() #
update() #
if first_click == 0: #
first_click = x #
else: #
sleep(1) #
draw(first_click,x) #
first_click = 0 #
wn = Screen()
draw(0,4)
wn.onclick(click)
wn.mainloop()
done()To find the attractors for each r, the code calculates the sequence elements until a repeating element is observed. Continuing the sequence from that point until another element is repeated ensures that all the attractors are visited. The code includes a tolerance mechanism to check if a number has already been observed. It allows a difference of 0.000001. This value can be decreased to obtain quicker results with less precision.
Here is the formed bifurcation diagram of the logistic map:
The program allows the user to click twice on the diagram at different points to focus on a narrower window of the diagram. Zooming in reveals strange arches and internal structures inside the chaotic region.
Within the chaos there are islands of order that appear as tiny repetitions of the initial structure of the diagram.
Summary
The mathematical expression xn+1=rxn(1−xn) is simple at first glance but surprising. Depending on r, it contains clearly defined regions as well as chaotic regions. With Python, one can explore the chaotic regions and discover surprising recurring patterns within them. The mathematical expression serves as a beautiful and easy-to-understand example of the concepts underlying chaos theory and fractal structures.