This is one of the 100 recipes of the IPython Cookbook, the definitive guide to high-performance scientific computing and data science in Python.

12.1. Plotting the bifurcation diagram of a chaotic dynamical system

  1. We import NumPy and matplotlib.
In [ ]:
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
  1. We define the logistic function by:
$$f_r(x) = rx(1-x)$$

Our discrete dynamical system is defined by the recursive application of the logistic function:

$$x_{n+1}^{(r)} = f_r(x_n^{(r)}) = rx_n^{(r)}(1-x_n^{(r)})$$
In [ ]:
def logistic(r, x):
    return r*x*(1-x)
  1. We will simulate this system for 10000 values of $r$ linearly spaced between 2.5 and 4. Of course, we vectorize the simulation with NumPy.
In [ ]:
n = 10000
r = np.linspace(2.5, 4.0, n)
  1. We will simulate 1000 iterations of the logistic map, and we will keep the last 100 iterations to display the bifurcation diagram.
In [ ]:
iterations = 1000
last = 100
  1. We initialize our system with the same initial condition $x_0 = 10^{-5}$.
In [ ]:
x = 1e-5 * np.ones(n)
  1. We will also compute an approximation of the Lyapunov exponent, for every value of $r$. The Lyapunov exponent is defined by:
$$\lambda(r) = \lim_{n \to \infty} \frac{1}{n} \sum_{i=0}^{n-1} \log\left| \frac{df_r}{dx}\left(x_i^{(r)}\right) \right|$$
In [ ]:
lyapunov = np.zeros(n)
  1. Now, we simulate the system and we plot the bifurcation diagram. The simulation only involves the iterative evaluation of the function $f$ on our vector $x$. Then, to display the bifurcation diagram, we draw one pixel per point $x_n^{(r)}$ during the last 100 iterations.
In [ ]:
for i in range(iterations):
    x = logistic(r, x)
    # We compute the partial sum of the Lyapunov exponent.
    lyapunov += np.log(abs(r-2*r*x))
    # We display the bifurcation diagram.
    if i >= (iterations - last):
        plt.plot(r, x, ',k', alpha=.04)
plt.xlim(2.5, 4);
plt.title("Bifurcation diagram");

# We display the Lyapunov exponent.
plt.plot(r[lyapunov<0], lyapunov[lyapunov<0] / iterations,
         ',k', alpha=.2);
plt.plot(r[lyapunov>=0], lyapunov[lyapunov>=0] / iterations,
         ',r', alpha=.5);
plt.xlim(2.5, 4);
plt.ylim(-2, 1);
plt.title("Lyapunov exponent");

The bifurcation diagram brings out the existence of a fixed point for $r<3$, then two and four equilibria... until a chaotic behavior when $r$ belongs to certain areas of the parameter space. We observe an important property of the Lyapunov exponent: it is positive when the system is chaotic (in red here).

You'll find all the explanations, figures, references, and much more in the book (to be released later this summer).

IPython Cookbook, by Cyrille Rossant, Packt Publishing, 2014 (500 pages).