{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> This is one of the 100 recipes of the [IPython Cookbook](http://ipython-books.github.io/), the definitive guide to high-performance scientific computing and data science in Python.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 12.1. Plotting the bifurcation diagram of a chaotic dynamical system"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "1. We import NumPy and matplotlib."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "2. We define the logistic function by:\n",
    "\n",
    "$$f_r(x) = rx(1-x)$$\n",
    "\n",
    "Our discrete dynamical system is defined by the recursive application of the logistic function:\n",
    "\n",
    "$$x_{n+1}^{(r)} = f_r(x_n^{(r)}) = rx_n^{(r)}(1-x_n^{(r)})$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def logistic(r, x):\n",
    "    return r*x*(1-x)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "3. 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."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "n = 10000\n",
    "r = np.linspace(2.5, 4.0, n)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "4. We will simulate 1000 iterations of the logistic map, and we will keep the last 100 iterations to display the bifurcation diagram."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "iterations = 1000\n",
    "last = 100"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "5. We initialize our system with the same initial condition $x_0 = 10^{-5}$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "x = 1e-5 * np.ones(n)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "6. We will also compute an approximation of the Lyapunov exponent, for every value of $r$. The Lyapunov exponent is defined by:\n",
    "\n",
    "$$\\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|$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "lyapunov = np.zeros(n)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "7. 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."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "plt.figure(figsize=(6,7));\n",
    "plt.subplot(211);\n",
    "for i in range(iterations):\n",
    "    x = logistic(r, x)\n",
    "    # We compute the partial sum of the Lyapunov exponent.\n",
    "    lyapunov += np.log(abs(r-2*r*x))\n",
    "    # We display the bifurcation diagram.\n",
    "    if i >= (iterations - last):\n",
    "        plt.plot(r, x, ',k', alpha=.04)\n",
    "plt.xlim(2.5, 4);\n",
    "plt.title(\"Bifurcation diagram\");\n",
    "\n",
    "# We display the Lyapunov exponent.\n",
    "plt.subplot(212);\n",
    "plt.plot(r[lyapunov<0], lyapunov[lyapunov<0] / iterations,\n",
    "         ',k', alpha=.2);\n",
    "plt.plot(r[lyapunov>=0], lyapunov[lyapunov>=0] / iterations,\n",
    "         ',r', alpha=.5);\n",
    "plt.xlim(2.5, 4);\n",
    "plt.ylim(-2, 1);\n",
    "plt.title(\"Lyapunov exponent\");\n",
    "plt.tight_layout();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "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.\n",
    "We observe an important property of the Lyapunov exponent: it is positive when the system is chaotic (in red here)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> You'll find all the explanations, figures, references, and much more in the book (to be released later this summer).\n",
    "\n",
    "> [IPython Cookbook](http://ipython-books.github.io/), by [Cyrille Rossant](http://cyrille.rossant.net), Packt Publishing, 2014 (500 pages)."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.4.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}