M62, TD4, Exo 2 : simulation

M67_TD4_exo2.md

Vous pouvez ouvrir ce gist dans le cloud :

• avec  : [OVH] [MyBinder] [Gesis] [Google Cloud] ;
• avec  : [Google Colab]

requirements.txt

numpy
matplotlib
scipy

TD4_exo2.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "hawaiian-switzerland",
   "metadata": {},
   "source": [
    "# TD4, exo 2 (simulation)\n",
    "\n",
    "## Chargement des bibliothèques standards"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "million-strength",
   "metadata": {},
   "outputs": [],
   "source": [
    "# numpy pour les calculs (vectoriels)\n",
    "import numpy as np\n",
    "# matplotlib (référencé comme `plt` ici) pour les graphiques\n",
    "from matplotlib import pyplot as plt\n",
    "# la bibliothèque pour tracer les trajectoires\n",
    "from scipy.integrate import odeint\n",
    "# la bibliothèque qui permet d'utiliser (le décorateur) @interact\n",
    "from ipywidgets import interact"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "saving-relief",
   "metadata": {},
   "source": [
    "## Une fonction utile\n",
    "\n",
    "La fonction `champ_normalise` permet de tracer le champs de vecteur $F(Y,t)$ **normalisé** associé au problème de Cauchy $Y'=F(Y,t)$ avec $Y : t \\mapsto (x(t),y(t)) \\in \\mathbb{R}$.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "aerial-market",
   "metadata": {},
   "outputs": [],
   "source": [
    "# affiche le champ normalisé sur la figure courante\n",
    "def champ_normalise(F, xmin, xmax, ymin, ymax, t=0, N=15):\n",
    "    X = np.linspace(xmin, xmax, N)  # abscisses des points de la grille\n",
    "    Y = np.linspace(ymin, ymax, N)  # ordonnées des points de la grille\n",
    "    U, V = F(np.meshgrid(X, Y), t)  # les composantes du champ de vecteurs\n",
    "    M = np.hypot(U, V)  # calcule la norme du vecteur (U,V)\n",
    "    M[M == 0] = 1  # évite la division par 0\n",
    "    U /= M  # normalise la composante U\n",
    "    V /= M  # normalise la composante V\n",
    "    return plt.quiver(X, Y, U, V, angles='xy')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "proper-poison",
   "metadata": {},
   "source": [
    "## La fonction F\n",
    "La fonction dans cet exercice dépend des paramètres $a,b,c,d$: \n",
    "$$ \n",
    "    F_{a,b,c,d}(Y,t) = (ax-bxy,-cy+dxy) \\text{ pour } Y=(x,y).\n",
    "$$\n",
    "\n",
    "Pour cette raison nous définisson la fonction `FF(a,b,c,d)` qui retourne une fonction `F(Y,t)` qui opère sur le 2-vecteur `Y`. Ainsi `Y[0]` représente la quantité $x$ et `Y[1]` représente la quantité $y$. \n",
    "\n",
    "*Remarque: Comme il s'agit d'une équation autonome ici, le paramètre `t` n'est pas utilisé pour le calcul de `F`.*"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "aerial-gambling",
   "metadata": {},
   "outputs": [],
   "source": [
    "# FF(a,b,c,d) retourne la fonction F\n",
    "FF = lambda a, b, c, d: lambda Y, t: [a * Y[0] - b * Y[0] * Y[1], -c * Y[1] + d * Y[0] * Y[1]]\n",
    "\n",
    "# test\n",
    "F = FF(1, 2, 3, 4)  # la fonction F en fcontion de r=1 et a=2\n",
    "F([3, 4], 0)  # la valeur de F en (3,4) à l'instant t=0"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fallen-funds",
   "metadata": {},
   "source": [
    "## Simulation\n",
    "\n",
    "On peut choisir les paramètres $a,b,c,d$, ainsi que la condition initiale $Y(0) = (x,y)$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "republican-exclusion",
   "metadata": {},
   "outputs": [],
   "source": [
    "# choix de la fenêtre\n",
    "xmin, xmax, ymin, ymax = -1, 10, -1, 10\n",
    "\n",
    "# discréisation du temps (100 instants dans [0,3])\n",
    "t = np.linspace(0, 3, 100)\n",
    "\n",
    "@interact(a=(0, 10), b=(0, 10), c=(0, 10), d=(0, 10), x=(0, xmax, .5), y=(0, xmax, .5))\n",
    "def dyn_soution(a=8, b=2, c=8, d=4, x=2, y=4):\n",
    "    # La fonction F en fonction des paramètres a,b,c,d\n",
    "    F = FF(a, b, c, d)\n",
    "\n",
    "    # on détermine la fenêtre d'affichage\n",
    "    plt.axis([xmin, xmax, ymin, ymax])\n",
    "    # on trace les isoclines (qui contiennent les axes)\n",
    "    plt.plot([xmin, xmax], [0, 0], \"r\")\n",
    "    if d != 0:\n",
    "        plt.plot([c / d, c / d], [ymin, ymax], \"r\")\n",
    "    plt.plot([0, 0], [ymin, ymax], \"m\")\n",
    "    if b != 0:\n",
    "        plt.plot([xmin, xmax], [a / b, a / b], \"m\")\n",
    "\n",
    "    # on dessine le champs normalisé\n",
    "    champ_normalise(F, xmin, xmax, ymin, ymax, t)\n",
    "\n",
    "    # la partie de la courbe pour t dans [0,3] (le futur)\n",
    "    v = odeint(F, [x, y], t)\n",
    "    plt.plot(v[:, 0], v[:, 1], \"b\")\n",
    "\n",
    "    # le point de départ (la condition initiale)\n",
    "    plt.plot(x, y, 'or')\n",
    "    # le titre\n",
    "    plt.title(f\"Pour a,b,c,d={a},{b},{c},{d} avec x(0)={x:4.1f}, y(0)={y:4.1f}\")\n",
    "    # et on affiche tout ça\n",
    "    plt.show()"
   ]
  }
 ],
 "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.7.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}