M62_TP3_NOM_Prenom.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": false,
    "editable": false
   },
   "source": [
    "![Département de Mathématiques](https://ktzanev.github.io/logolabopp/math-ulille/math-ulille_100.gif)\n",
    "\n",
    "# TP 3 - Les équations différentielles autonomes de dimension 2\n",
    "\n",
    "Dans ce troisième TP nous allons étudier numériquement des équations autonomes du premier ordre de dimension 2.<br>\n",
    "Il s'agit d'exemples vus dans la fiche B de TD.\n",
    "\n",
    "**Attention :** Pour que le code entre les exercices soit au maximum similaire et permettre de faire des copier/coller, nous utilisons plusieurs fois les mêmes noms pour des fonctions différentes. **Ceci est une très mauvaise pratique qui mène à des résultats imprévisibles si on n'exécute pas les cellules dans l'ordre.**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": false,
    "editable": false
   },
   "source": [
    "On commence par charger les bibliothèques « standards » dont on aura besoin."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "deletable": false,
    "editable": false
   },
   "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",
   "metadata": {
    "deletable": false,
    "editable": false
   },
   "source": [
    "## Quelques fonctions utiles"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": false,
    "editable": false
   },
   "source": [
    "### Le champ de vecteurs\n",
    "\n",
    "De façon similaire à la feuille de TP1, nous allons créer une fonction `champ_vectoriel`, mais avec les différences suivantes :\n",
    "- la fonction `F` passé en paramètre est la fonction correspondante à l'équation autonome de Cauchy $Y'=F(Y)$.\n",
    "- elle a un paramètre supplémentaire `normalise` qui indique s'il faut normaliser ou pas le champ de vecteurs.\n",
    "\n",
    "La fonction `F` passé en paramètre doit opérer sur un **vecteur** `Y` (de dimension $2\\times n$) qui représente les $n$ points de la grille où on affiche le champ de vecteurs."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "deletable": false,
    "editable": false
   },
   "outputs": [],
   "source": [
    "# affiche le champ de vecteurs sur la figure courante\n",
    "# - F : représente la fonction de Y'=F(Y)\n",
    "# - xmin, Xmax, ymin, ymax : détermine la fenêtre d'affichage\n",
    "# - normalise : si `False` le champ n'est pas normalisé\n",
    "# - step : est le pas de discrétisations en chaque direction\n",
    "def champ_vectoriel(F, xmin, xmax, ymin, ymax, normalise=True, step=.5, **kw):\n",
    "    X = np.arange(xmin, xmax, step)  # abscisses des points de la grille\n",
    "    Y = np.arange(ymin, ymax, step)  # ordonnées des points de la grille\n",
    "    U, V = F(np.meshgrid(X, Y))  # les composantes du champ de vecteurs\n",
    "    if normalise:\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', **kw)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": false,
    "editable": false
   },
   "source": [
    "#### Exemple\n",
    "\n",
    "Considérons le champ de vecteur de la fonction $F(x,y) = (x^2-y^2,2x)$. Le code suivant montre l'influence des paramètres `normalise` et `step` sur l'affichage du champ de vecteurs.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "deletable": false,
    "editable": false
   },
   "outputs": [],
   "source": [
    "F = lambda Y: np.array([Y[0] * Y[0] - Y[1] * Y[1], 2 * Y[0]])\n",
    "\n",
    "# choix de la fenêtre\n",
    "xmin, xmax, ymin, ymax = -4, 4, -2, 2\n",
    "\n",
    "# choix des paramètres et de la condition initiale grâce à des curseurs\n",
    "@interact(step=(.1, 1, .1), normalise=True)\n",
    "def dyn_exemple1(step, normalise):\n",
    "\n",
    "    # taille (maximale) de la figure en pouces (ici 14\" x 7\")\n",
    "    plt.figure(figsize=(14, 7))\n",
    "    # et on préserve les proportions des axes\n",
    "    plt.axis('scaled')\n",
    "    # on détermine la fenêtre d'affichage\n",
    "    plt.axis([xmin, xmax, ymin, ymax])\n",
    "\n",
    "    # l'abscisse\n",
    "    plt.plot([xmin, xmax], [0, 0], \"y\")\n",
    "    # l'ordonnée\n",
    "    plt.plot([0, 0], [ymin, ymax], \"y\")\n",
    "\n",
    "    # on dessine le champ (normalisé ou pas)\n",
    "    champ_vectoriel(F, xmin, xmax, ymin, ymax, normalise=normalise, step=step, width=.001)\n",
    "\n",
    "    # le titre\n",
    "    plt.title(\"Exemple de champ de vecteur\")\n",
    "    # et on affiche tout ça\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": false,
    "editable": false
   },
   "source": [
    "### Les courbes de niveau\n",
    "\n",
    "La fonction `lignes_niveaux` permet de tracer les courbes de niveaux de la fonction `h` passée en premier paramètre (où $h:\\mathbb{R}^2\\to\\mathbb{R}$). Nous allons utiliser cette fonction pour dessiner les courbes de niveaux d'un intégral premier ou d'une fonction de Liapounov, par exemple."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "deletable": false,
    "editable": false
   },
   "outputs": [],
   "source": [
    "# la fonction qui trace les lignes de niveaux de h qui passe par les points xy\n",
    "# - h : une fonction de R^2 dans R\n",
    "# - xmin, xmax, ymin, ymax : le limites du tracé\n",
    "# - xy : une liste de points de la forme [[x1,x2,...],[y1,y2,...]] par lesquelles passent les niveaux\n",
    "# - N : le nombre de points de la discrétisation de x et de y\n",
    "# - fmt : le format de la valeur de la courbe de niveau\n",
    "# - alpha : le niveau de transparence entre 0 (transparent) et 1 (opaque)\n",
    "def lignes_niveaux(h, xmin, xmax, ymin, ymax, xy, N=100, fmt='%.1f', alpha=.5):\n",
    "    C = np.unique(h(np.array(xy)))  # les valeurs de h aux points xy dans l'ordre croissant\n",
    "    x = np.linspace(xmin, xmax, N)\n",
    "    y = np.linspace(ymin, ymax, N)\n",
    "    z = h(np.meshgrid(x, y))  # les valeurs de h aux points (x,y)\n",
    "    lignes_niveau = plt.contour(x, y, z, levels=C, alpha=alpha)  # les courbes de niveaux\n",
    "    plt.clabel(lignes_niveau, C, fmt=fmt)  # on ajuste le format des étiquettes des courbes"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": false,
    "editable": false
   },
   "source": [
    "#### Exemple\n",
    "\n",
    "Considérons la fonction à deux paramètres $h_{a,b}(x,y) = ax^2+by^2$. Cet exemple illustre l'affichage des courbes de niveaux (qui sont des ellipses dans ce cas) qui passent par les points $(1,0)$, $(1,1)$ et $(0,1)$, ainsi que l'influence des paramètres `N` et `alpha`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "deletable": false,
    "editable": false
   },
   "outputs": [],
   "source": [
    "# La fonction hh(a,b) retourne la fonction h_{a,b}\n",
    "hh = lambda a, b: lambda Y: a * Y[0] * Y[0] + b * Y[1] * Y[1]\n",
    "\n",
    "# choix de la fenêtre\n",
    "xmin, xmax, ymin, ymax = -2, 2, -2, 2\n",
    "\n",
    "# choix des paramètres et de la condition initiale grâce à des curseurs\n",
    "@interact(a=(.1, 1, .1), b=(.1, 1, .1), N=(3, 100, 10), alpha=(0, 1, .1))\n",
    "def dyn_exemple2(a, b, N, alpha=1):\n",
    "\n",
    "    # taille (maximale) de la figure en pouces (ici 14\" x 7\")\n",
    "    plt.figure(figsize=(14, 7))\n",
    "    # et on préserve les proportions des axes\n",
    "    plt.axis('scaled')\n",
    "    # on détermine la fenêtre d'affichage\n",
    "    plt.axis([xmin, xmax, ymin, ymax])\n",
    "\n",
    "    # l'abscisse\n",
    "    plt.plot([xmin, xmax], [0, 0], \"y\")\n",
    "    # l'ordonnée\n",
    "    plt.plot([0, 0], [ymin, ymax], \"y\")\n",
    "\n",
    "    # la fonction h qui dépend des paramètres a et b\n",
    "    h = hh(a, b)\n",
    "\n",
    "    # les points par lesquels passent les courbes\n",
    "    xy = [[1, 1, 0], [0, 1, 1]]\n",
    "\n",
    "    # on dessine les courbes de niveaux\n",
    "    lignes_niveaux(h, xmin, xmax, ymin, ymax, xy, N=N, alpha=alpha)\n",
    "\n",
    "    # on dessine les points\n",
    "    plt.plot(xy[0], xy[1], \"or\")\n",
    "\n",
    "    # le titre\n",
    "    plt.title(\"Exemple de champ de vecteur\")\n",
    "    # et on affiche tout ça\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": false,
    "editable": false
   },
   "source": [
    "### Ajustements de F"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": false,
    "editable": false
   },
   "source": [
    "#### Rajout du paramètre temps\n",
    "\n",
    "Comme la fonction `odeint` exige deux paramètres, un vecteur et un paramètre temps `t`, on doit rajouter un `t` à la fonction `F` correspondant à l'équation autonome $Y'=F(Y)$. La fonction utilitaire `avec_t` fait ça : on lui passe une fonction à un paramètre `Y` et elle retourne une fonction à deux paramètres : `Y` et `t`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "deletable": false,
    "editable": false
   },
   "outputs": [],
   "source": [
    "avec_t = lambda F: lambda Y, t: F(Y)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": false,
    "editable": false
   },
   "source": [
    "### Ralentir le temps en dehors de la fenêtre d'affichage\n",
    "\n",
    "Si nous voulons afficher les trajectoires dans une fenêtre à priori nous ne savons pas dans quel intervalle il faut faire varier le temps. Si on donne un temps trop court, on n'affiche qu'un bout de l'orbite, si on donne un temps trop long la méthode `odeint` peut « exploser » (soit pour des raisons « théoriques » car il y a une explosion en temps fini, soit pour des raisons « numériques »).\n",
    "\n",
    "Pour éviter cela nous pouvons modifier `F` en dehors de la fenêtre d'affichage pour éviter ces problèmes, tout en gardant les même orbites dans la fenêtre d'affichage (où `F` reste non modifiée). \n",
    "\n",
    "C'est ce qu'on fait avec la fonction `ralentir(F, xmin, xmax, ymin, ymax)` qui retourne une version « ralentie » de `F` en dehors de la fenêtre `[xmin, xmax, ymin, ymax]`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "deletable": false,
    "editable": false
   },
   "outputs": [],
   "source": [
    "# C'est une fonction C^1 qui :\n",
    "#  - est égale à 1 dans le disque de rayon r\n",
    "#  - décroit exponentiellement en dehors de ce disque\n",
    "def cutoff(Y, r):\n",
    "    nY = np.maximum(Y[0] * Y[0] + Y[1] * Y[1] - r * r, 0)\n",
    "    return np.exp(-nY)\n",
    "\n",
    "# la fonction qui « ralentie » F en dehors de la fenêtre [xmin, xmax, ymin, ymax]\n",
    "def ralentir(F, xmin, xmax, ymin, ymax):\n",
    "    r = max(abs(xmax), abs(xmin), abs(ymax), abs(ymin)) * np.sqrt(2)  # le rayon d'un disque qui contient la fenêtre\n",
    "    return lambda Y: F(Y) * cutoff(Y, r)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": false,
    "editable": false
   },
   "source": [
    "## Exercice 4 (exemple)\n",
    "\n",
    "### 1) La fonction de Cauchy\n",
    "\n",
    "La fonction dans cet exercice ne dépend pas de paramètres :\n",
    "$$ \n",
    "    F(x,y) = (x^2-y, x)\n",
    "$$\n",
    "\n",
    "nous définissons cette fonction $F(Y)$ pour $Y=(x,y)$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "deletable": false,
    "editable": false
   },
   "outputs": [],
   "source": [
    "F = lambda Y: np.array([Y[0] * Y[0] - Y[1], Y[0]])\n",
    "\n",
    "# vérification\n",
    "assert np.all(F(np.array([[0, 1], [0, 2]])) == np.array([[0, -1], [0, 1]])), \"F(0,0)=(0,0) et F(1,2)=(-1,1)\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": false,
    "editable": false
   },
   "source": [
    "### 2) Le graphique « dynamique » avec les solutions approchées \n",
    "\n",
    "Ici, comme en TP1, nous allons afficher le champ de vecteurs et les solutions approchées obtenues par `odeint`. Comme ces solutions « explosent » numériquement dans certaines conditions, nous allons utiliser par défaut une version « ralentie » de `F`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "deletable": false,
    "editable": false
   },
   "outputs": [],
   "source": [
    "# choix de la fenêtre\n",
    "xmin, xmax, ymin, ymax = -4, 4, -1.5, 2\n",
    "\n",
    "# choix des paramètres et de la condition initiale grâce à des curseurs\n",
    "@interact(x_0=(xmin, xmax, .1), y_0=(ymin, ymax, .1), normalise=True, ralentie=True)\n",
    "# le dessin pour la condition initiale (x_0,y_0)\n",
    "def dyn_solution(x_0=0, y_0=0, normalise=True, ralentie=True):\n",
    "\n",
    "    # taille (maximale) de la figure en pouces (ici 14\" x 7\")\n",
    "    plt.figure(figsize=(14, 7))\n",
    "    # et on préserve les proportions des axes\n",
    "    plt.axis('scaled')\n",
    "    # on détermine la fenêtre d'affichage\n",
    "    plt.axis([xmin, xmax, ymin, ymax])\n",
    "\n",
    "    # l'abscisse\n",
    "    plt.plot([xmin, xmax], [0, 0], \"y\")\n",
    "    # l'isocline 0 (= l'ordonnée)\n",
    "    plt.plot([0, 0], [ymin, ymax], \"m\")\n",
    "    # l'isocline infinie (la parabole y = x^2)\n",
    "    u = np.linspace(xmin, xmax, 100)\n",
    "    plt.plot(u, u * u, \"r\")\n",
    "\n",
    "    # on dessine le champ (normalisé ou pas)\n",
    "    champ_vectoriel(F, xmin, xmax, ymin, ymax, normalise=normalise, width=.001)\n",
    "\n",
    "    # la solution particulière d'orbite la parabole y = x^2 - 1/2\n",
    "    u = np.linspace(2 * xmin, 2 * xmax, 100)\n",
    "    plt.plot(u / 2, u * u / 4 - 1 / 2, \"g\")\n",
    "\n",
    "    # ralentir ou pas\n",
    "    Fr = ralentir(F, xmin, xmax, ymin, ymax) if ralentie else F\n",
    "    # la fonction F avec le paramètre t en plus pour `odeint`\n",
    "    Ft = avec_t(Fr)\n",
    "    # le temps (future, le passé étant -t)\n",
    "    t = np.linspace(0, 5, 100)\n",
    "    # la partie de la courbe pour t dans [0,5] (le futur)\n",
    "    x, y = odeint(Ft, [x_0, y_0], t).transpose()\n",
    "    plt.plot(x, y, \"b\")\n",
    "    # la partie de la courbe pour t dans [-5,0] (le passé)\n",
    "    x, y = odeint(Ft, [x_0, y_0], -t).transpose()\n",
    "    plt.plot(x, y, \"g--\")\n",
    "\n",
    "    # le point de départ (la condition initiale)\n",
    "    plt.plot(x_0, y_0, 'or')\n",
    "    # le titre\n",
    "    plt.title(f\"La solution qui vérifie x(0)={x_0:4.1f}, y(0)={y_0:4.1f}\")\n",
    "    # et on affiche tout ça\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": false,
    "editable": false
   },
   "source": [
    "## Exercices 1 et 7 (exemple)\n",
    "\n",
    "### 1) La fonction de l'équation de Cauchy\n",
    "\n",
    "Dans la deuxième partie de l'exercice 1 de la fiche B nous avons étudié un système de deux équations deux inconnues dépend de deux paramètres $a > 0$ et $r > 0$ qui se réduit au problème de Cauchy autonome :\n",
    "$$\n",
    "    F_{r,a}(Y) = (-rSI,rSI-aI) \\text{ pour } Y=(S,I).\n",
    "$$\n",
    "L'exercice 7 étudie le même système avec $a=1$ et $r=2$, mais il rajoute l'étude de une intégrale première et d'une fonction de Liapounov.\n",
    "\n",
    "Pour cette raison nous définissons la fonction `FF(r, a)` qui retourne une fonction `F(Y)` qui opère sur le 2-vecteur `Y`. Ainsi `Y[0]` représente la quantité $S$ et `Y[1]` représente la quantité $I$ dans l'exercice 1, ou $x$ et $y$ respectivement de l'exercice 7."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "deletable": false,
    "editable": false
   },
   "outputs": [],
   "source": [
    "# FF(a,b) retourne la fonction F\n",
    "# On peut le faire avec `lambda` ainsi :\n",
    "# FF = lambda r, a: lambda Y: np.array([-r * Y[0] * Y[1], r * Y[0] * Y[1] - a * Y[1]])\n",
    "# Mais on peut le faire aussi avec `def` de façon plus « éloquente »\n",
    "def FF(r, a):\n",
    "    def F(Y):\n",
    "        S, I = Y  # S = Y[0], I = Y[1]\n",
    "        Fx = -r * S * I\n",
    "        Fy = r * S * I - a * I\n",
    "        return np.array([Fx, Fy])\n",
    "\n",
    "    return F\n",
    "\n",
    "# test\n",
    "assert np.all(FF(1, 2)([3, 4]) == np.array([-12, 4])), \"Pour r=1 et a=2, la fonction F vaut [-12, 4] en [3, 4].\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": false,
    "editable": false
   },
   "source": [
    "### 2) Une intégrale première\n",
    "\n",
    "Une intégrale première pour $S>0$ (et $I>0$) est\n",
    "$$\n",
    "    h_{a,r}(S,I) = I + S - \\frac{a}{r}\\log(S).\n",
    "$$\n",
    "*Vérifier qu'il s'agit bien d'un intégral premier.*\n",
    "\n",
    "**Remarque :** Le cas limite de Lotka-Volterra (exercice 2) avec $a=0$ et $b=d$ coïncide avec le système $S$ et $I$ étudié ici. Ainsi l'intégrale première peut être déduite de la question 7 de ce même exercice. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "deletable": false,
    "editable": false
   },
   "outputs": [],
   "source": [
    "# la fonction hh construit et retourne la fonction h sur la base des paramètres a,b,c,d.\n",
    "def hh(r, a):\n",
    "    def h(Y):\n",
    "        S, I = Y\n",
    "        return I + S - np.log(S) * a / r\n",
    "\n",
    "    return h"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": false,
    "editable": false
   },
   "source": [
    "Pour représenter les courbes de niveau de la fonction `h = hh(a,r)` nous allons utiliser la fonction définie plus haut `lignes_niveaux` comme suit (avec $a=1$ et $r=2$, comme dans l'exercice 7)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "deletable": false,
    "editable": false
   },
   "outputs": [],
   "source": [
    "# choix de la fenêtre\n",
    "xmin, xmax, ymin, ymax = -.1, 3, -.1, 3\n",
    "\n",
    "# choix des paramètres du problème\n",
    "r = 2\n",
    "a = 1\n",
    "\n",
    "# taille (maximale) de la figure en pouces (ici 14\" x 7\")\n",
    "plt.figure(figsize=(14, 7))\n",
    "# et on préserve les proportions des axes\n",
    "plt.axis('scaled')\n",
    "# on détermine la fenêtre d'affichage\n",
    "plt.axis([xmin, xmax, ymin, ymax])\n",
    "# on trace les axes (qui coïncident avec l'isocline infinie)\n",
    "plt.plot([xmin, xmax], [0, 0], \"r\")  # l'abscisse = les points fixes\n",
    "plt.plot([0, 0], [ymin, ymax], \"r\")\n",
    "\n",
    "# l'isocline I_0 (sans les points fixes = abscisse)\n",
    "plt.plot([a / r, a / r], [ymin, ymax], \"m\")\n",
    "# les courbes de niveaux\n",
    "h = hh(r, a)\n",
    "yl = np.linspace(0, ymax, 11)[1:]  # 10 nombres entre 0 (exclu) et ymax\n",
    "xl = np.ones_like(yl) * a / r  # x = a/r\n",
    "e = 1e-5  # presque 0, pour éviter la valeur 0 dans le logarithme de h\n",
    "lignes_niveaux(h, e, xmax, 0, ymax, [xl, yl])\n",
    "\n",
    "# on dessine le champ normalisé\n",
    "F = FF(r, a)\n",
    "champ_vectoriel(F, 0, xmax, 0, ymax, step=.25, width=.003)\n",
    "\n",
    "# le titre\n",
    "plt.title(f\"Le champ de vecteur et l'intégrale première pour a={a}, r={r}\")\n",
    "# et on affiche tout ça\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": false,
    "editable": false
   },
   "source": [
    "### 3) Une fonction de Liapounov\n",
    "\n",
    "La fonction de Liapounov de l'exercice 7 (pour $a=1$ et $r=2$, mais en réalité elle ne dépend pas de ces paramètres) est la fonction :\n",
    "$$\n",
    "    l(S, I) = S + I.\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "deletable": false,
    "editable": false
   },
   "outputs": [],
   "source": [
    "liapunov = lambda Y: Y[0] + Y[1]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": false,
    "editable": false
   },
   "source": [
    "Pour représenter les courbes de niveaux de cette fonction de Liapounov on peut procéder comme avec l'intégrale première."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": false,
    "editable": false
   },
   "source": [
    "### 4) Le graphique « dynamique » avec les solutions approchées\n",
    "\n",
    "Nous somme prêt maintenant d'intégrer tout ça dans une image « dynamique » qui contient :\n",
    "- le champ de vecteurs (normalisé ou pas);\n",
    "- les isoclines $0$ et $\\infty$ que nous avons trouvé en TD ;\n",
    "- le point de départ $(t_0, y_0)$ ;\n",
    "- l'orbite de la solution approchée obtenue avec `odeint` qui vérifie $y(t_0) = y_0$.\n",
    "\n",
    "**Remarque :** Dans cet exemple nous n'avons pas d'explosion numérique (les orbites sont bornées), donc nous n'avons pas besoin de « ralentir » `F`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "deletable": false,
    "editable": false
   },
   "outputs": [],
   "source": [
    "# choix de la fenêtre\n",
    "xmin, xmax, ymin, ymax = -.1, 10, -.1, 3\n",
    "\n",
    "# options possibles\n",
    "options_champ = {\"ne pas afficher\": None, \"normalisé\": True, \"non normalisé\": False}\n",
    "options_niveaux = {\"ne pas afficher\": None, \"fonction de Liapounov\": 1, \"intégrale première\": 2}\n",
    "\n",
    "# choix des paramètres et de la condition initiale grâce à des curseurs\n",
    "@interact(r=(1, 10), a=(1, 10), x_0=(0, xmax, .5), y_0=(0, ymax, .5), champ=options_champ, niveaux=options_niveaux)\n",
    "# le dessin pour la condition initiale (x_0,y_0)\n",
    "def dyn_solution(r=1, a=3, x_0=4, y_0=2, champ=True, niveaux=2):\n",
    "\n",
    "    # taille (maximale) de la figure en pouces (ici 14\" x 7\")\n",
    "    plt.figure(figsize=(14, 7))\n",
    "    # et on préserve les proportions des axes\n",
    "    plt.axis('scaled')\n",
    "    # on détermine la fenêtre d'affichage\n",
    "    plt.axis([xmin, xmax, ymin, ymax])\n",
    "\n",
    "    # on trace les axes (qui coïncident avec l'isocline infinie)\n",
    "    plt.plot([xmin, xmax], [0, 0], \"r\")  # l'abscisse = les points fixes\n",
    "    plt.plot([0, 0], [ymin, ymax], \"r\")\n",
    "    # l'isocline I_0 (sans les points fixes = abscisse)\n",
    "    plt.plot([a / r, a / r], [ymin, ymax], \"m\")\n",
    "\n",
    "    # les courbes de niveaux\n",
    "    if niveaux is not None:\n",
    "        # soit de l'inégale première, soit de Liapounov\n",
    "        g = liapunov if niveaux == 1 else hh(r, a)\n",
    "\n",
    "        # on dessine les niveaux qui passe par des points (a/r,yl)\n",
    "        yl = np.linspace(0, ymax, 10)[1:]\n",
    "        xl = np.ones_like(yl) * a / r\n",
    "        # pour éviter log(0) dans l'intégrale première nous commençons S à 10^-5\n",
    "        lignes_niveaux(g, 1e-5, xmax, 0, ymax, [xl, yl])\n",
    "\n",
    "    F = FF(r, a)\n",
    "    # on dessine (ou pas) le champ normalisé\n",
    "    if champ is not None:\n",
    "        champ_vectoriel(F, 0, xmax, 0, ymax, normalise=champ, step=.25, width=.001)\n",
    "\n",
    "    # la fonction F avec le paramètre t en plus pour `odeint`\n",
    "    Ft = avec_t(F)\n",
    "    # le temps\n",
    "    t = np.linspace(0, 3, 100)\n",
    "    # la partie de la courbe pour t dans [0,3] (le futur)\n",
    "    x, y = odeint(Ft, [x_0, y_0], t).transpose()\n",
    "    plt.plot(x, y, \"b\")\n",
    "    # la partie de la courbe pour t dans [-3,0] (le passé)\n",
    "    x, y = odeint(Ft, [x_0, y_0], -t).transpose()\n",
    "    plt.plot(x, y, \"g--\")\n",
    "\n",
    "    # le point de départ (la condition initiale)\n",
    "    plt.plot(x_0, y_0, 'or')\n",
    "    # le titre\n",
    "    plt.title(f\"Pour a={a}, r={r}, la solution qui vérifie avec S(0)={x_0:4.1f}, I(0)={y_0:4.1f}\")\n",
    "    # et on affiche tout ça\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": false,
    "editable": false
   },
   "source": [
    "## Exercice 2 (à faire)\n",
    "\n",
    "Dans cet exercice nous avons étudié le système différentiel dit de Lotka-Volterra, où a, b, c, d paramètres\n",
    "(strictement) positifs :\n",
    "$$\n",
    "  \\begin{align}\n",
    "    \\dot{x} &= \\phantom{-}ax - bxy\\\\\n",
    "    \\dot{y} &= -cy + dxy\n",
    "  \\end{align}\n",
    "$$\n",
    "\n",
    "\n",
    "### 1) La fonction de Cauchy\n",
    "\n",
    "Définissez une fonction `FF(a, b, c, d)` qui en fonction des paramètres `a,b,c,d` retourne la fonction `F` correspondante à l'équation autonome de Cauchy $Y'=F(Y)$ avec $Y=(x,y)$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "deletable": false
   },
   "outputs": [],
   "source": [
    "# La fonction FF(a,b,c,d) retourne la fonction F qui prend un 2-vecteur comme paramètre\n",
    "\n",
    "# test\n",
    "F = FF(1, 2, 3, 4)  # la fonction F en fonction des paramètres a=1, b=2, c=3, d=4\n",
    "Y = np.array([[0, 1], [0, 2]])  # vérification pour les vecteurs (0,0) et (1,2)\n",
    "R = np.array([[0, -3], [0, 2]])  # le résultat attendu de F(Y)\n",
    "assert np.all(F(Y) == R), \"Il faut avoir F(0,0)=(0,0) et F(1,2)=(-3,2) pour ces paramètres\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": false,
    "editable": false
   },
   "source": [
    "### 2) Intégrale première\n",
    "\n",
    "Définissez une fonction `hh(a, b, c, d)` qui en fonction des paramètres `a,b,c,d` retourne la fonction `h` correspondante à l'intégrale première trouvée en TD :\n",
    "$$\n",
    "    h(x,y) = a\\log(y)-by+c\\log(x)-dx.\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "deletable": false
   },
   "outputs": [],
   "source": [
    "# La fonction hh(a,b,c,d) retourne la fonction h qui prend un 2-vecteur Y=(x,y) comme paramètres\n",
    "\n",
    "# test\n",
    "h = hh(0, 1, 2, 3)  # la fonction h pour les paramètres a=0, b=1, c=2, d=3\n",
    "Y = np.ones((2, 3))  # Y = [[1,1,1],[1,1,1]]\n",
    "x, y = Y\n",
    "assert np.all(h(Y) == -x - 3 * y), \"Il faut avoir F(1,1)=-4 pour ces paramètres\"\n",
    "Y = np.array([2, 3])\n",
    "assert f\"{hh(1, 2, 3, 4)(Y):.2f}\" == \"-10.82\", f\"{hh(1, 2, 3, 4)(Y):.2f} == -10.82 ?\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": false,
    "editable": false
   },
   "source": [
    "### 3) Le graphique « dynamique » avec les solutions approchées\n",
    "\n",
    "Inspirer vous des exemples traités pour afficher un graphique dynamique qui : \n",
    "- permet de choisir les paramètres `a,b,c,d` entre 0 et 10 ;\n",
    "- permet de choisir la condition initiale `x_0, y_0` entre 0 et 10 ;\n",
    "- qui affiche plusieurs courbes de niveaux de l'intégrale première `h = hh(a, b, c, d)`;\n",
    "- qui (quand c'est possible) affiche les isoclines $I_0$ et $I_\\infty$ ;\n",
    "- qui affiche le champ de vecteurs normalisé ;\n",
    "- contient la solution approchée obtenue avec `odeint`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "deletable": false
   },
   "outputs": [],
   "source": [
    "# choix de la fenêtre\n",
    "xmin, xmax, ymin, ymax = -1, 10, -1, 10\n",
    "\n",
    "# choix des paramètres et de la condition initiale grâce à des curseurs\n",
    "@interact(a=(0, 10), b=(0, 10), c=(0, 10), d=(0, 10), x_0=(0, xmax, .5), y_0=(0, xmax, .5))\n",
    "# le dessin avec la solution pour la condition initiale (x_0,y_0)\n",
    "def dyn_solution(a=7, b=2, c=7, d=2, x_0=3.5, y_0=3.5):\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": false,
    "editable": false
   },
   "source": [
    "### 4) Observations\n",
    "\n",
    "1. Que se passe t-il quand $a=0$ ? Reconnaissez-vous un système étudié ?\n",
    "1. Et quand $a=0$ et $c=0$ qu'observe t-on ? Comment on peut expliquer ça avec la signification « biologique » des paramètres ? "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": false
   },
   "source": [
    "**Solution :**\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": false,
    "editable": false
   },
   "source": [
    "## Exercice 3 (à faire)\n",
    "\n",
    "### 1) La fonction de Cauchy\n",
    "\n",
    "La fonction du problème de Cauchy étudié dans cet exercice ne dépend pas de paramètres :\n",
    "$$ \n",
    "    F(Y) = \\big(x(1+\\frac{3}{x^2+y^2+2}),-y(1-\\frac{3}{x^2+y^2+2})\\big) \\text{ pour } Y=(x,y).\n",
    "$$\n",
    "\n",
    "Définissez cette fonction."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "deletable": false
   },
   "outputs": [],
   "source": [
    "# Définition de la fonction F\n",
    "\n",
    "# test\n",
    "assert np.all(F([[0, 1], [1, 0]]) == np.array([[0, 2], [0, 0]])), \"F(0,1) = (0,0) et F(1,0) = (2,0)\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": false,
    "editable": false
   },
   "source": [
    "### 2) Le graphique « dynamique » avec les solutions approchées\n",
    "\n",
    "Inspirer vous des exemples traités pour afficher un graphique dynamique qui : \n",
    "- permet de choisir la condition initiale `x_0, y_0` entre 0 et 10 ;\n",
    "- affiche les isoclines $I_0$ et $I_\\infty$ (trouvées en TD);\n",
    "- affiche le champ de vecteurs (normalisé ou pas en fonction du paramètre `normalise`)\n",
    "- contient la solution approchée obtenue avec `odeint`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "deletable": false
   },
   "outputs": [],
   "source": [
    "# choix de la fenêtre\n",
    "xmin, xmax, ymin, ymax = -.1, 11, -.1, 4\n",
    "\n",
    "# choix de la condition initiale et du paramètre de normalisation\n",
    "@interact(x_0=(0, xmax, .25), y_0=(0, ymax, .25), normalise=True)\n",
    "# le dessin avec la solution pour la condition initiale (x_0,y_0)\n",
    "def dyn_solution(x_0=0, y_0=0, normalise=True):\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": false,
    "editable": false
   },
   "source": [
    "### 3) Observations\n",
    "\n",
    "De combien d'orbite est composé la partie positive de l'ordonnée d'après les observations ? (décrivez votre expérience) "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": false
   },
   "source": [
    "**Solution :**\n"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}