Period-doubling bifurcations leading to chaos, and their bifurcation diagrams.sagews

Note: If you want to use any of the interactives in this worksheet (highly recommended), click the link that says “Open in CoCalc” at the top right of the page.

The period-doubling path to chaos

and the corresponding bifurcation diagrams

We have seen in class that in the discrete-time logistic model, $$X_{t+1} = r \cdot X_t \cdot (1 - X_t) ,$$ as we vary the parameter $r$ between $0$ and $4$, we get a series of bifurcations.

When $0 \le r < 3$, the system goes to a stable equilibrium. That equilibrium point is at $0$ when $r \le 1$, but is nonzero when $1 < r < 3$. Right at $r = 3$, a Hopf bifurcation occurs, so that when $r > 3$ (up to about $3.45$) the system goes to a stable oscillation between two $X$ values.

The interactive below allows you to explore this.

@interact
def logistic(r=slider(0, RDF(3.43), RDF(0.01), default=2, label="$r$")):
    f(X) = r*X*(1 - X)
    solution = [0.1]
    for n in xrange(50):
        solution.append(f(solution[-1]))
    p = list_plot(solution, pointsize=30)
    p += list_plot(solution, plotjoined=True)
    p.show(ymin=0, ymax=1, axes_labels=("$t$", "$X$"))

Interact: please open in CoCalc

For $r$ values a little larger than $3.45$, the system starts to display a new type of behavior: after some transient behavior, it settles into a cycle that repeats itself, but rather than cycling between just two values of $X$, it cycles between four values! That is, instead of just repeating a pattern like A-B-A-B-A-B... (where A represents a low value of $X$ and B represents a higher value), it repeats a pattern of the form A-C-B-D-A-C-B-D-A-C-B-D. This is still periodic behavior, because it is repeating the same values over and over, but now the period is $4$ rather than $2$, because it takes $4$ time steps before it repeats.

Since increasing the parameter $r$ past the threshold of $3.45$ caused the oscillations to go from period $2$ to period $4$, we say that a period-doubling bifurcation occurred at $r \approx 3.45$.

This new “period $4$ behavior” can be seen in the steady-state part of the graph below.

r = 3.5
f(X) = r*X*(1 - X)
solution = [0.1]
for n in xrange(50):
    solution.append(f(solution[-1]))
p = list_plot(solution, pointsize=30)
p += list_plot(solution, plotjoined=True)
p += text("$r = {}$".format(RDF(r)), (25, 1), color="black", fontsize=16)
p.show(ymin=0, ymax=1, axes_labels=("$t$", "$X$"))

Slightly farther along, after about $r = 3.54$, the pattern of behavior changes again, this time to A-E-D-G-B-F-C-H-A-E-D-G-B-F-C-H-A-E-D-G-B-F-C-H. That is, the cycle now has period $8$. Another period-doubling bifurcation has occurred!

After around $r = 3.56$, the behavior changes to period $16$ (another period-doubling bifurcation), and after about $r = 3.565$, it changes to period $32$ (yet another period-doubling bifurcation). You can explore all of these in the interactive below.

What you should be noticing about the $r$ values mentioned above is that they are getting closer and closer together. Here is a summary:

• $r = 3.00$: Hopf bifurcation, start of period 2 oscillations

• $r \approx 3.45$: First period-doubling bifurcation, start of period 4 oscillations

• $r \approx 3.543$: Second period-doubling bifurcation, start of period 8 oscillations

• $r \approx 3.562$: Third period-doubling bifurcation, start of period 16 oscillations

• $r \approx 3.5665$: Fourth period-doubling bifurcation, start of period 32 oscillations

• etc...

In fact, between the Hopf bifurcation at $r = 3$ and somewhere around $r = 3.56995$, an infinite number of period-doubling bifurcations occur! But because they keep occuring closer and closer together, all infinitely many of them fit within the relatively short interval of $r$ values from $3.45$ to $3.56995$.

Then, for values of $r$ greater than about $3.56995$, the behavior changes again, this time to chaotic behavior.

In the interactive below, you can see the behavior for specific $r$ values in this range, plus a few larger $r$ values showing chaotic behavior. The first plot shows the time series from $t = 200$ to $t = 250$. The first $200$ iterations have been discarded, to get rid of any transient behavior and just show the steady-state pattern. Although it's a discrete-time model, this plot was made with plotjoined=True to help show the pattern of the oscillation. In the same interactive, the second plot also shows only the steady-state behavior (by discarding the first $200$ iterations again), but this time shows many more iterations, with plotjoined=False and very small points. This makes it much easier to see how many distinct $X$ values are hit in the oscillation pattern, so it's easy to see what the period of the oscillation is.

@interact
def logistic(r=slider([3.4, 3.5, 3.56, 3.566, 3.575, 3.6, 4])):
    f(X) = r*X*(1 - X)
    solution = [0.1]
    for t in xrange(1, 2000):
        solution.append(f(solution[-1]))
    p = list_plot(zip(range(200, 250), solution[200:250]), plotjoined=True)
    p.show(axes_labels=("$t$", "$X$"), ymin=0, ymax=1)
    p = list_plot(zip(range(200, 2000), solution[200:2000]), pointsize=1)
    p.show(axes_labels=("$t$", "$X$"), ymin=0, ymax=1)

Interact: please open in CoCalc

We can create a bifurcation diagram that summarizes all of these behaviors in a way that's very similar to the bifurcation diagrams we studied in LS 30A. The parameter ($r$) value goes on the horizontal axis, and the vertical axis is the state variable ($X$), so that each vertical slice of the graph represents a copy of the state space for a specific $r$ value. However, instead of just plotting the equilibrium points (and indicating their stability with colors or a line pattern), here we will plot, for each $r$ value, all of the $X$ values that are part of the steady-state oscillation pattern at that value of $r$. So counting the points in any vertical slice of the diagram will tell you the period of the oscillation. Here is that diagram:

As you can see, the Hopf bifurcation at $r = 3$ shows up here as a “branching”, where the curve forks into two curves. Likewise, each period-doubling bifurcation appears as an $r$ value at which all of the curves simultaneously branch in the same way: $2$ curves branch into $4$ curves, then $4$ curves branch into $8$ curves, then $8$ curves branch into $16$ curves, etc.

The regions of chaotic behavior (most of the graph where $r$ is greater than about $3.6$) show up as a sort of fuzzy gray area, because for these $r$ values, the system no longer generates a finite sequence of $X$ values that repeats periodically. Instead, since chaotic behavior is aperiodic, the sequence of states will include infinitely many different $X$ values. To produce this picture, instead of just plotting a point for every $X$ value, we have shaded each point lighter or darker depending on how frequently the sequence of $X$ values hit that point. When the behavior is periodic, every $X$ value in the periodic pattern gets hit just as often as every other one, so those points are colored uniformly black. But when it's chaotic, some $X$ values occur more frequently than others, so some are shaded darker and some are lighter.

(Note: I have not included the code used to produce these diagrams in this worksheet, but they were produced in CoCalc using SageMath. If you'd like to see/use the code yourself, contact me.)

Here is a higher-resolution version, showing just the more interesting $r$ range from $2.5$ to $4$.

One thing that you might notice in these diagrams is that, for some $r$ values greater than $3.6$, the chaotic behavior seems to suddenly disappear for a little while, and a few individual curves appear again, indicating that the system has returned to periodic behavior! Probably the most visible of these regions is around $r = 3.83$. What does the period seem to be there? If you look closely, shortly after that $r$ value, the system goes through another period-doubling bifurcation. In fact, although it's hard to see it without zooming in much more, the system will once again undergo infinitely many more period-doubling bifurcations (all between $r = 3.83$ and about $r = 3.85$) before returning to chaotic behavior! This means that if you were to zoom in on a rectangle defined by roughly $3.83 \le r \le 3.85$ and $0.9 \le X \le 1$, the diagram would look very similar to the region of the graph from $r = 2.5$ to $r = 3.6$... that is, most of the whole graph that you see above! Mathematicians call this property self-similarity, and images that are self-similar are often referred to as fractals. The self-similarity of bifurcation diagrams like this one is one of the deep connections between chaotic dynamical systems and fractal geometry.

The region we've just described isn't the only one with periodic behavior in the middle of all the chaos. Everywhere in the chaotic region where you see a white vertical band, there should be periodic behavior with some finite period. With the help of the diagram above, see if you can use the interactive below to find $r$ values where the behavior is periodic with period $3$, period $6$, period $5$, and period $7$. What other periods can you find?

@interact
def logistic(r=slider(RDF(3.569), 4, RDF(0.001), label="$r$")):
    f(X) = r*X*(1 - X)
    solution = [0.1]
    for t in xrange(1, 2000):
        solution.append(f(solution[-1]))
    p = list_plot(zip(range(200, 250), solution[200:250]), plotjoined=True)
    p.show(axes_labels=("$t$", "$X$"), ymin=0, ymax=1)
    p = list_plot(zip(range(200, 2000), solution[200:2000]), pointsize=1)
    p.show(axes_labels=("$t$", "$X$"), ymin=0, ymax=1)

Interact: please open in CoCalc

What about other chaotic systems?

Does a similar pattern of bifurcations happen in other models that show chaotic behavior? And in particular, can we find anything equivalent to the period-doubling bifurcations that we saw above in a continuous-time (differential equation) model? Could we even make a bifurcation diagram in a similar way, and if so, what will it look like?

The answers to all of these questions are yes, yes, yes, and surprisingly similar to the bifurcation diagrams that we just saw!

My favorite example of this in continuous time is the Rössler model (which you covered in one of your labs as the “Romeo-Juliet-Tybalt” model). As a reminder, in this model, $R$ represents Romeo's love (or hatred if negative) for Juliet, $J$ represents Juliet's love/hatred for Romeo, and $T$ represents Tybalt's hatred of Romeo. The system of differential equations is $$\begin{cases} R' &= -J - T \\ J' &= R + 0.1 J \\ T' &= 0.1 - a T + R T \end{cases}$$ In your lab, you studied this model with the parameter $a$ (in the $T'$ equation) set to $14$. But since we're interested in studying bifurcations of this model, we need to vary a parameter, so the parameter $a$ is the one we'll work with.

Note that since there are three state variables in this model (as is necessary to get chaotic behavior from an ODE), the graphs of trajectories shown below will all be 3D graphs.

The first graph below shows, in different colors, a few different trajectories of this model for the parameter value $a = 2$. You should see a very familiar feature here.

a = 2
vectorfield(R, J, T) = (-J - T, R + 0.1*J, 0.1 - a*T + R*T)
t_range = srange(0, 100, 0.01)
p = axes3d((R, -5, 5), (J, -5, 5), (T, 0, 5), grid=False)
solution = desolve_odeint(vectorfield, (0, 4, 0), t_range, [R, J, T])
p += list_plot(solution, plotjoined=True, color="red")
solution = desolve_odeint(vectorfield, (1, 1, 0), t_range, [R, J, T])
p += list_plot(solution, plotjoined=True, color="fuchsia")
solution = desolve_odeint(vectorfield, (0, -3.78, 0.034), t_range, [R, J, T])
p += list_plot(solution, plotjoined=True, color="blue", thickness=2)
p.show(aspect_ratio=(1, 1, 1), frame=False)

3D rendering not yet implemented

As you can see, the slightly thicker blue loop, which the red and pink trajectories approach asymptotically, is a limit cycle attractor. So when $a = 2$, this model already shows periodic behavior. (It turns out that it's not possible to get a stable equilibrium point for any value of the parameter $a$, so we can't get a Hopf bifurcation by varying this parameter.)

From now on, since we are only interested in the steady-state behavior, we will do as we did for the discrete-time logistic model above, and throw away the earlier part of our simulation to get rid of any transients.

The next two graphs show just the limit cycle attractor for the parameter values $a = 4$ and $a = 6$. What happens between these two values?

t_range = srange(0, 100, 0.01)
a = 4
vectorfield1(R, J, T) = (-J - T, R + 0.1*J, 0.1 - a*T + R*T)
solution1 = desolve_odeint(vectorfield1, (1, 1, 0), t_range, [R, J, T])
a = 6
vectorfield2(R, J, T) = (-J - T, R + 0.1*J, 0.1 - a*T + R*T)
solution2 = desolve_odeint(vectorfield2, (1, 1, 0), t_range, [R, J, T])

start = int(75 / 0.01)
show("$\Large a = 4$")
p = axes3d((R, -12, 12), (J, -12, 12), (T, 0, 5), grid=False)
p += list_plot(solution1[start:], plotjoined=True, color="blue")
p.show(aspect_ratio=(1, 1, 1), frame=False)

print("\n\n\n")

show("$\Large a = 6$")
p = axes3d((R, -12, 12), (J, -12, 12), (T, 0, 5), grid=False)
p += list_plot(solution2[start:], plotjoined=True, color="blue")
p.show(aspect_ratio=(1, 1, 1), frame=False)

$\Large a = 4$

3D rendering not yet implemented

$\Large a = 6$

3D rendering not yet implemented

What just happened?

If you look closely, each of the above trajectories is a closed loop, so the behavior we're seeing here is periodic. However, in the second one, the trajectory has wound around the $T$ axis twice before closing the loop. If you consider the amount of time required for the trajectory to wind around the $T$ axis once, it's not hard to imagine (even though we can't see $t$ in the graph) that it would take about twice as much time to go around twice. In other words, it's likely the the period of the oscillatory behavior in the second ($a = 6$) graph is about twice as long as the period of the oscillation in the first ($a = 4$) graph.

To summarize what we just said: as the parameter $a$ was increased, somewhere between $a = 4$ and $a = 6$ the period of the oscillation doubled. What does that mean? You guessed it. That's a period-doubling bifurcation!

To see what we mean about the period as a time interval, take a look at the two time series below for the $T$ variable. The period of the oscillation can be measured here as the horizontal distance between the peaks. In the first ($a = 4$) graph, that period is approximately $6$. In the second ($a = 6$) graph, it's approximately $12$. Hence, the period doubled.

p = list_plot(zip(t_range, solution1[:,2])[start:], plotjoined=True)
p += text("$a = 4$", (87.5, 9), color="black", fontsize=16)
p.show(ymax=10, axes_labels=("$t$", "$T$"))

p = list_plot(zip(t_range, solution2[:,2])[start:], plotjoined=True)
p += text("$a = 6$", (87.5, 9), color="black", fontsize=16)
p.show(ymax=10, axes_labels=("$t$", "$T$"))

The interactive below will show you both the trajectory and the time series (for the $T$ variable) of this model, for a handful of specific values of $a$. These values of $a$ show a sequence of period-doubling bifurcations, and eventually chaotic behavior.

@interact
def rossler(a=slider([4, 6, 8.5, 8.7, 8.73, 8.74, 8.75, 8.76, 9], label="$a$")):
    vectorfield(R, J, T) = (-J - T, R + 0.1*J, 0.1 - a*T + R*T)
    t_range = srange(0, 1000 if a < 8.74 else 1500, 0.01)
    start = int(900/0.01)
    solution = desolve_odeint(vectorfield, (1, 1, 0), t_range, [R, J, T])
    p = axes3d((R, -20, 20), (J, -20, 20), (T, 0, 20), grid=False)
    p += list_plot(solution[start:], plotjoined=True, color="blue")
    p.show(aspect_ratio=(1, 1, 1), frame=False)
    p2 = list_plot(zip(t_range, solution[:,2])[start:], plotjoined=True)
    p2.show(ymax=20, axes_labels=("$t$", "$T$"))

Interact: please open in CoCalc

As with the discrete-time logistic model that we considered earlier, we can create a bifurcation diagram showing the changes in behavior as we vary the parameter $a$. Here it is:

Notice how similar this is to the bifurcation diagram from the discrete-time logistic model. Of course, the precise details are different. But the overall shape or pattern of behavior is exactly the same: we see a sequence of period-doubling bifurcations, occurring closer and closer together, until suddenly chaotic behavior happens. Then, within the region of chaotic behavior, there are small windows of periodic behavior that appear for short parameter intervals, and even within these more period-doubling bifurcations can be seen. Also, within the region of chaotic behavior, some values of the state variable ($J$ in this case) occur more frequently, and others less frequently, resulting in darker or lighter shading. From this, a complex pattern can be seen through the chaos.

If you're following closely how these bifurcation diagrams are created, you might be wondering how this one can be made, since there are now $3$ state variables in this model instead of $1$, and the trajectory moves smoothly through the state space instead of jumping from one point to another as in a discrete-time model. In this case, we have “discretized” the steady-state behavior of the system by looking at the points where the trajectory intersects the $JT$ plane. We have then decreased these $3$-dimensional points to $1$-dimensional data by looking only at the value of the $J$ variable at each point (since the $T$ values at these points are quite close to $0$ and don't vary much). This is the reason that the vertical axis in the bifurcation diagram above is labeled $J$: it gives the $J$ coordinates of the points at which the trajectory intersects the $JT$ plane. Similar tricks can be employed in other continuous-time models with chaotic behavior.

You can use the interactive below to explore this model further for the whole range of values of $a$ shown above. With the help of the above bifurcation diagram, see if you can use this interactive to find $a$ values where the limit cycle attractor winds around the $T$ axis $3$ times before forming a loop. What about $6$ times, or $5$ times, or $7$ times?

@interact
def rossler(a=slider(1, 18, RDF(0.01), default=6, label="$a$")):
    vectorfield(R, J, T) = (-J - T, R + 0.1*J, 0.1 - a*T + R*T)
    t_range = srange(0, 1500, 0.01)
    start = int(900/0.01)
    solution = desolve_odeint(vectorfield, (1, 1, 0), t_range, [R, J, T])
    p = axes3d((R, -25, 25), (J, -25, 25), (T, 0, 40), grid=False)
    p += list_plot(solution[start:], plotjoined=True, color="blue")
    p.show(aspect_ratio=(1, 1, 1), frame=False)
    p2 = list_plot(zip(t_range, solution[:,2])[start:], plotjoined=True)
    p2.show(ymax=40, axes_labels=("$t$", "$T$"))

Interact: please open in CoCalc

The Hastings $3$-species food chain model

Recall the following model from class (and your textbook), in which $Z$ is the population of a top predator such as wolves, $Y$ is the population of a middle species, such as sheep, that is preyed upon by the top predator but is also a predator of the bottom species, $X$, which might be, for example, grass. $$\begin{cases} X' = r \cdot X \cdot \big( 1 - \frac{X}{k} \big) - m_1 \cdot \frac{X}{h_1 + X} \cdot Y \\ Y' = b_1 \cdot \big( m_1 \cdot \frac{X}{h_1 + X} \cdot Y \big) - d_1 \cdot Y - m_2 \cdot \frac{Y}{h_2 + Y} \cdot Z \\ Z' = b_2 \cdot \big( m_2 \cdot \frac{Y}{h_2 + Y} \cdot Z \big) - d_2 \cdot Z \end{cases}$$

We have seen in class that the parameter values $r = 1$, $k = 1$, $m_1 = 1.6$, $h_1 = 0.3$, $b_1 = 1$, $d_1 = 0.4$, $m_2 = 0.05$, $h_2 = 0.5$, $b_2 = 1$, $d_2 = 0.01$ gives rise to chaotic behavior. We will now look at what happens as we vary the parameter $m_1$, leaving all the other parameters the same.

Below is an interactive that allows you to explore this, followed by the resulting bifurcation diagram. Again, see if you can find values of the parameter that give rise to an ordinary limit cycle attractor, and to a limit cycle attractor that winds around $2$, $3$, or $4$ times before forming a closed loop.

@interact
def hastings(m1=slider(0.65, 1.8, RDF(0.01), label="$m_1$")):
    vectorfield(X, Y, Z) = (X*(1 - X) - m1*X*Y/(0.3 + X), m1*X*Y/(0.3 + X) - 0.4*Y - 0.05*Y*Z/(0.5 + Y), 0.05*Y*Z/(0.5 + Y) - 0.01*Z)
    initial_state = (12.5*sqrt(0.002704 - 0.0008*m1) + 0.4, 0.1, 12.5*m1 + 375*sqrt(0.002704 - 0.0008*m1) - 24.4)
    t_range = srange(0, 5000, 0.1)
    start = int(1000/0.1)
    solution = desolve_odeint(vectorfield, initial_state, t_range, [X, Y, Z])
    p = axes3d((X, 0, 1), (Y, 0, 0.5), (Z, 0, 11), grid=False, origin=(0, 0, 0))
    p += list_plot(solution[start:], plotjoined=True)
    p.show(aspect_ratio=(1, 2, 0.25), frame=False)

Interact: please open in CoCalc

Bifurcation diagram for the Hastings model (Note: The parameter is called $b_1$ in the diagram, but it corresponds to what we've called $m_1$ above.)

It is also possible to generate a bifurcation diagram like this for the Lorenz model (which was covered in your lab as “Romeo, Juliet, and Juliet's Nurse”). Once again, it will have most of the exact same features as all of the bifurcation diagrams we have looked at above. The precise shape will be different, but in terms of its overall look, you could say all of these look “about the same”.

Summary: What these bifurcation diagrams mean

In brief, a bifurcation diagram like the one above means that as we increase the parameter $r$, the steady-state behavior of the system changes in the following way:

$$\text{equilibrium} \to \text{oscillations} \to \text{complex oscillations} \to \text{chaos}$$

Here, “complex oscillations” refer to the oscillatory behavior with higher periods, after one or more period-doubling bifurcations have occurred.

Recall that the change from equilibrium to oscillations is a Hopf bifurcation, and the change from regular oscillations to complex oscillations is a period-doubling bifurcation. Also, within the region of complex oscillations, infinitely many more period-doubling bifurcations occur, closer and closer together. This is sometimes referred to as a “period-doubling cascade”, which then leads to chaotic behavior. In general, the above pattern is often referred to as the “period-doubling route to chaos”.