Interactive Lorentz System.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "## Exploring the Lorentz System\n",
    "\n",
    "The following code is [available on the web](http://ipywidgets.readthedocs.io/en/latest/examples/Lorenz%20Differential%20Equations.html)\n",
    "\n",
    "Another great visualization is found [here](https://jakevdp.github.io/blog/2013/02/16/animating-the-lorentz-system-in-3d/)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "from ipywidgets import interact, interactive\n",
    "from IPython.display import clear_output, display, HTML\n",
    "import numpy as np\n",
    "from scipy import integrate\n",
    "\n",
    "from matplotlib import pyplot as plt\n",
    "from mpl_toolkits.mplot3d import Axes3D\n",
    "from matplotlib.colors import cnames\n",
    "from matplotlib import animation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "def solve_lorenz(N=10, angle=0.0, max_time=4.0, sigma=10.0, beta=8./3, rho=28.0):\n",
    "\n",
    "    fig = plt.figure()\n",
    "    ax = fig.add_axes([0, 0, 1, 1], projection='3d')\n",
    "    ax.axis('off')\n",
    "\n",
    "    # prepare the axes limits\n",
    "    ax.set_xlim((-25, 25))\n",
    "    ax.set_ylim((-35, 35))\n",
    "    ax.set_zlim((5, 55))\n",
    "\n",
    "    def lorenz_deriv(x_y_z, t0, sigma=sigma, beta=beta, rho=rho):\n",
    "        \"\"\"Compute the time-derivative of a Lorenz system.\"\"\"\n",
    "        x, y, z = x_y_z\n",
    "        return [sigma * (y - x), x * (rho - z) - y, x * y - beta * z]\n",
    "\n",
    "    # Choose random starting points, uniformly distributed from -15 to 15\n",
    "    np.random.seed(1)\n",
    "    x0 = -15 + 30 * np.random.random((N, 3))\n",
    "\n",
    "    # Solve for the trajectories\n",
    "    t = np.linspace(0, max_time, int(250*max_time))\n",
    "    x_t = np.asarray([integrate.odeint(lorenz_deriv, x0i, t)\n",
    "                      for x0i in x0])\n",
    "\n",
    "    # choose a different color for each trajectory\n",
    "    colors = plt.cm.viridis(np.linspace(0, 1, N))\n",
    "\n",
    "    for i in range(N):\n",
    "        x, y, z = x_t[i,:,:].T\n",
    "        lines = ax.plot(x, y, z, '-', c=colors[i])\n",
    "        plt.setp(lines, linewidth=2)\n",
    "\n",
    "    ax.view_init(30, angle)\n",
    "    plt.show()\n",
    "\n",
    "    return t, x_t"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "#t, x_t = solve_lorenz(angle=0, N=10)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "613ddff78a4f451aa41181e90bd09155"
      }
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "w = interactive(solve_lorenz, angle=(0.,360.), max_time=(0.1, 4.0),\n",
    "                N=(0,50), sigma=(0.0,50.0), rho=(0.0,50.0))\n",
    "display(w)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": true,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "t, x_t = w.result"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{'N': 21,\n",
       " 'angle': 109.8,\n",
       " 'beta': -0.06667,\n",
       " 'max_time': 2.3,\n",
       " 'rho': 23.3,\n",
       " 'sigma': 13.1}"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "w.kwargs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": true,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "xyz_avg = x_t.mean(axis=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(21, 3)"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "xyz_avg.shape"
   ]
  }
 ],
 "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.6.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}