Chaos and the Logistic Difference Equation

The logistic difference equation (or logistic map) is defined as follows:

This equation, and more complicated variations of it, have been studied by ecologists since the 1950s. It can be used as a simple model of the population dynamics of many animals like fish and insects. Given a population x_n at the nth year and a value for the growth parameter r, this equation returns x_n+1, the population for the next year.

What is especially interesting about the logistic difference equation though is the qualitative change in the behaviour of the system for different values of the parameter r. Depending on what we set r to be, our modelled population may go extinct. Or it might settle at a particular size and stay there forever. Or it could bounce between two sizes, or between four sizes, or even bounce all over the place chaotically!

So as a first post for the new year (and the new blog!), we'll use some of the gorgeous Racket plotting facilities to explore this logistic difference equation and its chaotic behaviour.

The Logistic Difference Equation

The logistic difference equation is given by x_n+1=rx_n(1-x_n). We can translate this straight into a Racket function like this:


As mentioned, the equation is used as a very simple model of an animal population that changes year by year. Here x_n is a value between 0 (extinction) and 1 (maximum possible population). So then the (1-x_n) term is also between 0 and 1, and works to limit the population growth – when x_n gets large, (1-x_n) gets small, and vice versa.

The parameter r represents the growth rate, chosen to be a value between 0 and 4.
(Why might this be so? Well, we want x_n=1 to represent the maximum possible population size. So if we look at the logistic equation as a continuous function, say y=rx(1-x), then calculus says that the maximum of this function occurs at dy/dx=r-2rx=0; that is, x=1/2, for a value of y(1/2)=r/4.

But we would want to use this value of y as the next value of x_n+1 in the case of the difference equation. So ensuring that r stays between 0 and 4 keeps us within the 0 ≤ x_n+1 ≤ 1 bounds.)

So with the function coded, we can explore its behaviour by playing around with it and seeing what happens to the population for different values of the growth rate parameter r.

Graphing Populations for Different Values of r

Now the logistic difference equation tells us exactly how to generate the next value for the population size given the current value. For a chosen value of the growth rate parameter r and a starting point x₀ in the interval (0,1), we just need to repeatedly apply the equation to generate a series of population sizes, which we can then graph against time.

So we can write a Racket function iterate-logistic-equation that takes as arguments the number of times to generate a new population value (the steps), the growth rate r, and x-init for the starting point. The function then iteratively applies the logistic equation to the most recent value, stores them all in a list and then constructs vectors for the graph points by pairing each population size with its time index.

The function then returns a list of these points (vectors), which we can pass straight to the line plotter in the Racket plot library. We can wrap all of this up in a tidy function plot-logistic-equation that takes the number of steps, the growth rate r and the starting point x-init and displays the graph for x(n). The code for all of this is shown below.


With the code written, we simply call it at the REPL to generate graphs of the logistic difference equation for different values of r:
 
r=0.1   (plot-logistic-equation 20 0.1 0.4)

r=0.6   (plot-logistic-equation 20 0.6 0.4)

r=1.1   (plot-logistic-equation 20 1.1 0.4)

r=2.1   (plot-logistic-equation 20 2.1 0.4)

r=2.8   (plot-logistic-equation 20 2.8 0.4)

As we can sort of see from the plots, for 0≤r<1, the population size drops to zero and stays there. That is, for values of the growth rate less than 1, we have one steady state: a population size of 0 (extinction.)

As r increases above 1 though, we still have a single steady state, but this time it is not zero; the population size converges to some nonzero value and stays in that equilibrium.

So far there's nothing too weird going on here. However, suppose we keep increasing the growth rate r...

r=2.9   (plot-logistic-equation 50 2.9 0.4)

r=3.0   (plot-logistic-equation 50 3.0 0.4)

r=3.2   (plot-logistic-equation 50 3.2 0.4)

Here we are starting to see oscillation between two different population sizes. We've had a bifurcation – for some value of the parameter r between 2.9 and 3.2, the qualitative behaviour of the system changes, and instead of a single steady state we now have two steady states that the population size alternates between. (So we've had a period-doubling bifurcation: now we have two time steps between repeated population sizes.)

We can increase r even further and see what happens... 

r=3.4   (plot-logistic-equation 50 3.4 0.4)

r=3.5   (plot-logistic-equation 50 3.5 0.4)

Now we've had another period-doubling bifurcation, and the population size oscillates between four different steady states.

But if we keep increasing the growth rate r now, we see something really strange... 

r=3.55  (plot-logistic-equation 60 3.55 0.4)

r=3.65  (plot-logistic-equation 60 3.65 0.4)

r=3.75  (plot-logistic-equation 60 3.75 0.4)

r=3.90  (plot-logistic-equation 60 3.90 0.4)

At this point, the population size appears to bounce around erratically, not settling on any specific values at all. We have chaos!

So what happened? We should absolutely investigate this phenomenon. And as it turns out, we have some quick and simple analytic tools that we can use to get a sense of what is going on.

Cobweb Diagrams for the Logistic Difference Equation

One thing that we can do to visualise how the system moves towards steady states (or doesn't, as the case may be) is to plot the cobweb diagram for the system. As detailed by Sternberg (2010), we plot a graph of the logistic difference equation y=rx(1-x) on the x-y plane, together with the diagonal line y=x. Then for our given starting point x₀, we do the following:

1. Draw the line connecting (x₀,0) and (x₀,y(x₀)). That is, we draw a line up from the x-axis at x₀ until we reach the graph of y=rx(1-x).
2. Our relation x_n+1=rx_n(1-x_n) specifies that this value of y(x₀) is actually the next value for x₁. So we draw a line horizontally across to the diagonal y=x. We'll be at the point (y(x₀),y(x₀)) now (ie. We've 'set' x₁=y(x₀).)
3. We now draw the vertical line connecting (x₁,x₁) to (x₁,y(x₁)), and continue on in this fashion to get each of the points x_n. As we'll see, in this way we draw the iterations for the population sizes.

We can code up this strategy in a Racket function plot-cobweb-diagram, which takes the number of steps, the growth rate parameter r and the starting point x-init and generates the set of vertical and horizontal lines for the cobweb diagram, plotting them together with the functions y=rx(1-x) and y=x on the x-y plane. The code for this is as follows:


We can plot the cobweb diagrams for different values of the growth parameter r and different starting points x-init in a straight-forward way from the Racket REPL.
  
r=0.5   (plot-cobweb-diagram 30 0.5 0.7)

r=0.9   (plot-cobweb-diagram 30 0.9 0.7)

In this way, we can easily see how the convergence to 0 could happen for 0≤r<1. We can also see how the single steady states might be reached:
 
r=1.6   (plot-cobweb-diagram 30 1.6 0.8)

r=2.4   (plot-cobweb-diagram 30 2.4 0.1)

 
Now recall that when we plotted the logistic difference equation values before, we seemed to have a bifurcation occurring somewhere between r=2.9, r=3.0 and r=3.2, and we started seeing oscillation between two different steady states. A look at the cobweb diagrams for these values of the growth parameter may help us see what's going on:
  
r=2.9   (plot-cobweb-diagram 200 2.9 0.4)

r=3.0   (plot-cobweb-diagram 200 3.0 0.4)

r=3.2   (plot-cobweb-diagram 200 3.2 0.4)

From the diagrams, we can see that the oscillation for the two steady states seems to coincide with the succession of vertical and horizontal lines in the cobweb diagram converging to a square. We can see a similar sort of thing happening for the four steady state cases:
  
r=3.4   (plot-cobweb-diagram 200 3.4 0.4)

r=3.5   (plot-cobweb-diagram 200 3.5 0.4)

And the cobweb diagrams also show us what happens with with the iterations when we cross that threshold into chaotic behaviour:

r=3.55  (plot-cobweb-diagram 200 3.55 0.4)

r=3.65  (plot-cobweb-diagram 200 3.65 0.4)

r=3.75  (plot-cobweb-diagram 200 3.75 0.4)

r=3.90  (plot-cobweb-diagram 200 3.90 0.4)

Here the vertical and horizontal lines don't converge to a point, or a square, or anything like that. Instead we've got a spiral that bounces all over the place!

Bifurcation Diagram for the Logistic Difference Equation

We saw in the previous plots and diagrams that the steady state population values changed as we increased the growth parameter r. So perhaps it might be enlightening to graph the steady state values against that parameter, so we can see the points at which the bifurcations occur?

Such a diagram is called a bifurcation diagram, and we'll generate one for the logistic difference equation with the following method, also detailed by Sternberg (2010):

• We'll slowly vary the bifurcation parameter r across the range [0,4].
• For each value of the parameter r, we'll iterate from a starting point for some number of steps (say 100), and generate a list of the population sizes we hit. Since we usually jump around for a bit at the start of the iteration, we'll keep track of only the last 30 or so points.
• We will then plot these steady state values against the value of the growth parameter r.

Note that this method will also allow us to see the sort of chaotic behaviour that happens for large values of the growth parameter, since we keep the last 30 point each time for the plot. These won't technically be 'steady states' obviously, but we'll plot them anyway.

We can write another Racket function plot-bifurcation-diagram to do all this, which takes an initial point x-init as its argument and displays the bifurcation diagram. The function looks like this:

 
We can then generate the bifurcation diagram by calling the function from the REPL as always. Here the value for x-init doesn't actually matter much (provided it's between 0 and 1 of course), so we'll arbitrarily choose x-init to be 0.3.

   (plot-bifurcation-diagram 0.3)

Here we can clearly see how the steady state is 0 for growth rates between 0 and 1, and that the steady state population size increases with r afterward, up until r is approximately 3.0. At that point we have our first period-doubling bifurcation: we can see how the curve splits into two separate branches.

We then have another period-doubling bifurcation at about r=3.5, where the two branches each split again. And as r keeps increasing, we quickly see the onset of chaos, represented by the scattered mess of points on the right-hand side.

Interestingly enough, we can notice from the bifurcation diagram that we appear to have a thin sliver of white in the otherwise blue chaos region at about r=3.85. Sure enough, when we plot the logistic difference equation and the cobweb diagram for that value of the growth rate, we do appear to have some slightly more predictable behaviour: 

r=3.85   (plot-logistic-equation 100 3.85 0.3)

r=3.85   (plot-cobweb-diagram 200 3.85 0.3)

So that was a quick (and hopefully not too intense!) introduction to the chaotic behaviour that can be found in the logistic difference equation, and perhaps also an introduction to the fascinating field of Chaos Theory. Until next time!

[For convenience, the complete Racket source file for all of the code in this blog post can be found here.]

References

Gleick, J. 1988, Chaos: Making a New Science, 1st edn, William Heinemann Ltd, London, UK.

Sternberg, S. 2010, Dynamical Systems, 1st edn, Dover Publications Inc, Mineola, NY.