03-scipy-integrate.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Modèles d'intégration avec la bibliothèque `scipy`\n",
    "\n",
    "La bibliothèque `scipy.integrate` propose un modèle d'intégrateur. Il suffit de fournir une fonction qui transforme un état en sa dérivée à l'instant $t$ pour que le module intègre cette fonction.\n",
    "\n",
    "Prenons l'exemple du lancer de boulet de canon. On considère une boule de masse $m$ et de rayon $r$ lancée du point $x_0$ avec une vitesse initiale $v_0$.\n",
    "\n",
    "On commence par négliger toutes les forces autres que celle exercée par la gravité. Le champ est uniforme avec $g=9,81\\,m\\,s^{-2}$.\n",
    "\n",
    "Le principe fondamental de la dynamique donne le système d'équations:\n",
    "\\begin{cases}\n",
    "\\;\\ddot{x}(t) =  -g\\cdot\\vec{z}\\\\\n",
    "\\;\\dot{x}(0) =  v_0\\\\\n",
    "\\;x(0) =  x_0\n",
    "\\end{cases}\n",
    "\n",
    "Ce système s'intègre simplement en un polynôme du second degré. On propose ici de résoudre ce système à l'aide de l'intégrateur fourni dans `scipy`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "\n",
    "from scipy.integrate import odeint\n",
    "\n",
    "def forces(state, t):\n",
    "    m = 20.                     # kilograms\n",
    "    g = np.array([0, -9.8])    # metres per second\n",
    "    \n",
    "    dstate = state.copy()\n",
    "    dstate[:2] = state[2:]\n",
    "    dstate[2:] = g\n",
    "        \n",
    "    return dstate\n",
    "\n",
    "# x_0, z_0, \\dot{x}_0, \\dot{z}_0\n",
    "state0 = np.array([0., 100., 100., 100.])\n",
    "\n",
    "t = np.arange(0.0, 30.0, 0.1)\n",
    "\n",
    "state = odeint(forces, state0, t)\n",
    "\n",
    "state=state[np.where(state[:,1] >=0)]\n",
    "\n",
    "plt.plot(state[:,0], state[:,1])\n",
    "\n",
    "plt.xlabel('x')\n",
    "plt.ylabel('z')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-warning\">\n",
    "**Exercice:** Ajouter la force de frottement à l'air et comparer les trajectoires sur le même graphe.\n",
    "</div>\n",
    "\n",
    "**Note**: La force de frottement à l'air s'exerce en opposition au vecteur vitesse. Son module s'exprime de la manière suivante:\n",
    "$$F = \\frac{1}{2}\\cdot C_x \\cdot \\rho \\cdot S \\cdot V^{2}$$\n",
    "\n",
    "On prendra $C_x = 0.5$ (sphère), une surface de référence $S =0,1\\,m^2$ et $\\rho = 1,184\\;kg\\,m^{-3}$ la masse volumique de l'air."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# %load solutions/ballistic.py\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### L'attracteur de Lorenz\n",
    "\n",
    "L'attracteur de Lorenz est une bonne illustration d'un système chaotique. Il est modélisé par les équations suivantes:\n",
    "\n",
    "$$\\dot{x} = \\sigma (y - x)$$\n",
    "$$\\dot{y} = (\\rho - z) x - y$$\n",
    "$$\\dot{z} = x\\, y - \\beta\\, z$$\n",
    "\n",
    "On prendra $\\sigma = 10$, $\\rho = 28$ et $\\beta = \\dfrac{8}{3}$\n",
    "\n",
    "<div class=\"alert alert-warning\">\n",
    "**Exercice:** Intégrer les équations sur l'intervalle donné.\n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.integrate import odeint\n",
    "\n",
    "σ = 10\n",
    "ρ = 28\n",
    "β = 8/3\n",
    "\n",
    "def lorenz(state, t):\n",
    "    x, y, z = state\n",
    "\n",
    "    # fill here\n",
    "    \n",
    "\n",
    "    return [1, 1, 1]\n",
    "\n",
    "state0 = [2.0, 3.0, 4.0]\n",
    "t = np.arange(0.0, 30.0, 0.01)\n",
    "state = odeint(lorenz, state0, t)\n",
    "\n",
    "# do some fancy 3D plotting\n",
    "from mpl_toolkits.mplot3d import Axes3D\n",
    "\n",
    "fig = plt.figure(figsize=(10, 7))\n",
    "ax = fig.gca(projection='3d')\n",
    "ax.plot(*state.T, alpha=0.4)\n",
    "ax.set_xlabel('x')\n",
    "ax.set_ylabel('y')\n",
    "ax.set_zlabel('z')\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Le système semble avoir deux équilibres stables : les deux circuits. Le système gravite autour d'un voisinage pour parfois changer de voisinage. Le nombre de cycles dans un voisinage et le moment auquel le système change de voisinage sont extrêmement dépendants des conditions initiales."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-warning\">\n",
    "**Exercice:** Afficher la trace précédent dans le plan 2D. Puis afficher les traces (position et trace pendant les 100 dernières millisecondes) de deux points extrêmement proches pour détecter le moment où les trajectoires divergent.\n",
    "</div>\n",
    "\n",
    "On pourra afficher les traces toutes les 5 secondes, puis toutes les secondes entre 15 et 20 secondes."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# %load solutions/lorenz.py\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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": 2
}