M62_TP2_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 2 - Les équations différentielles de dimension 1 (deuxième partie)\n",
    "\n",
    "Dans ce deuxième TP nous continuons l'étudier numériquement des équations du premier ordre de dimension 1.<br>\n",
    "Il s'agit d'un exemple vu dans les fiches A de TD, ainsi que deux exercices du DS2 de 2021.\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": [
    "## Une fonction utile\n",
    "\n",
    "On reprend du TP1 la fonction `champ_normalise` qui permet de tracer le champ de vecteur $(1,f(y,t))$ **normalisé** associé au problème de Cauchy $y'=f(y,t)$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "deletable": false,
    "editable": false
   },
   "outputs": [],
   "source": [
    "# affiche le champ normalisé sur la figure courante\n",
    "# - f représente la fonction de y'=f(y,t)\n",
    "# - tmin, tmax, ymin, ymax est la fenêtre d'affichage\n",
    "# - N est le nombre de discrétisations en chaque direction\n",
    "def champ_normalise(f, tmin, tmax, ymin, ymax, N=15, **kw):\n",
    "    F = lambda Y: [np.ones_like(Y[0]), f(Y[1], Y[0])]  # la fonction du problème autonome associé\n",
    "    T = np.linspace(tmin, tmax, 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(T, Y))  # 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(T, Y, U, V, angles='xy', **kw)  # trace le champ de vecteurs sur la grille de NxN points"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": false,
    "editable": false
   },
   "source": [
    "## Exercice 11 (à faire)\n",
    "\n",
    "\n",
    "Notre équation autonome étant\n",
    "$$\n",
    "    y'=\\frac{y}{2}\\big( 1-\\frac{3}{y^2+2}\\big)\n",
    "$$\n",
    "la fonction $f$ à définir est $f(y,t) =\\frac{y}{2}\\big( 1-\\frac{3}{y^2+2}\\big)$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "deletable": false
   },
   "outputs": [],
   "source": [
    "# la fonction f\n",
    "\n",
    "# test\n",
    "assert f(0, 0) == 0, \"f(0,*) = 0\"\n",
    "assert f(1, 1) == 0, \"f(1,*) = 0\"\n",
    "assert np.shape(f(np.ones(3), np.zeros(3))) == (3, ), \"la fonction doit accepter des vecteurs comme paramètres\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": false,
    "editable": false
   },
   "source": [
    "### La représentation graphique\n",
    "\n",
    "Nous avons vu en TD que ce problème possède trois solutions constantes : $1$, $0$ et $-1$. Et que entre ces constantes la solution est monotone, car la dérivée est de signe constant.\n",
    "\n",
    "Adapter les représentations graphiques des exercices du TP1 à cet exercice, représentant les solutions constantes ainsi que la solution qui vérifie $y(t_0) = y_0$.\n",
    "\n",
    "*Aide pour la question d'après :* Si possible, représenter sur l'ordonné l'orbite $y(\\mathbb{R})\\subset \\mathbb{R}$ décrite par cette solution."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "deletable": false
   },
   "outputs": [],
   "source": [
    "# choix de la fenêtre\n",
    "tmin, tmax, ymin, ymax = -10, 10, -5, 5\n",
    "\n",
    "# choix de la condition initiale grâce à des curseurs\n",
    "@interact(t_0=(tmin, tmax, .5), y_0=(ymin, ymax, .25))\n",
    "# le dessin pour la condition initiale (t_0,y_0)\n",
    "def dyn_solution(t_0=0, y_0=0):\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": false,
    "editable": false
   },
   "source": [
    "**Question :** En observant le graphique, l'ensemble $\\mathbb{R}$ est la réunion de combien d'orbites ?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": false
   },
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": false,
    "editable": false
   },
   "source": [
    "## Exercice 1 du DS2 de 2021 (à faire)\n",
    "\n",
    "On considère le problème de Cauchy, avec $|\\cdot|$ qui est la valeur absolue sur $\\mathbb{R}$ :\n",
    "\n",
    "$$\n",
    "    \\begin{cases}\n",
    "        y'(t)=|y(t)+t|\\\\\n",
    "        y(0)=0\n",
    "    \\end{cases}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": false,
    "editable": false
   },
   "source": [
    "1) Établir que ce problème admet une unique solution $y$ de classe $C^1$ globale, i.e. définie sur $\\mathbb{R}$ tout entier."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": false
   },
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": false,
    "editable": false
   },
   "source": [
    "2) Calculer cette solution."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": false
   },
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": false,
    "editable": false
   },
   "source": [
    "## La solution globale\n",
    "Définir la fonction `s` qui est la solution globale vérifiant `s(0)==0`. Cette fonction doit pouvoir prendre comme argument des `np.array`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "deletable": false
   },
   "outputs": [],
   "source": [
    "# Définition de la fonction s(t), solution qui vérifie la condition initiale y(0)=0\n",
    "\n",
    "# test\n",
    "assert s(0) == 0, \"La fonction doit vérifier la condition initiale s(0) == 0\"\n",
    "assert s(-2) < s(1), \"La fonction doit être strictement croissante\"\n",
    "assert np.shape(s(np.ones(3))) == (3, ), \"la fonction doit accepter des vecteurs comme paramètres\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": false,
    "editable": false
   },
   "source": [
    "## La représentation graphique\n",
    "\n",
    "Adapter les représentations graphiques dynamiques des exercices du TP1 à cet exercice, représentant la solution globale `s` ainsi que la solution trouvée par `odeint` qui vérifie $y(t_0) = y_0$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "deletable": false
   },
   "outputs": [],
   "source": [
    "# Définition de la fonction f(y,t) de l'équation y'=f(y,t)\n",
    "\n",
    "# test\n",
    "assert f(0, 0) == 0, \"f(0,0) = 0\"\n",
    "assert f(1, -1) == 0, \"f(1,-1) = 0\"\n",
    "assert np.shape(f(np.ones(3), np.zeros(3))) == (3, ), \"la fonction doit accepter des vecteurs comme paramètres\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "deletable": false
   },
   "outputs": [],
   "source": [
    "# choix de la fenêtre\n",
    "tmin, tmax, ymin, ymax = -2, 2, -2, 2\n",
    "\n",
    "# choix de la condition initiale grâce à des curseurs\n",
    "@interact(t_0=(tmin, tmax, .5), y_0=(ymin, ymax, .25))\n",
    "# le dessin pour la condition initiale (t_0,y_0)\n",
    "def dyn_solution(t_0=0, y_0=0):\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": false,
    "editable": false
   },
   "source": [
    "## Exercice 3 du DS2 de 2021 (à faire)\n",
    "\n",
    "On va terminer ce TP par un exercice dont les solutions maximales ne sont pas globales, car elles explosent en temps fini.\n",
    "\n",
    "Soit $y$ la solution maximale du problème de Cauchy :\n",
    "\n",
    "$$\n",
    "    \\begin{cases}\n",
    "        y'(t)=y(t)+y^3(t)\\\\\n",
    "        y(0)=1\n",
    "    \\end{cases}\n",
    "$$\n",
    "\n",
    "On peut voir facilement que ce problème vérifie les conditions de Cauchy-Lipschitz local, mais pas global. Et comme la fonction constante $0$  est une solution (globale) de cette équation autonome, on conclut que les autres solutions ne s'annulent pas.\n",
    "Puis en posant $z(t)= \\frac{1}{y^2(t)}$ et voit que $z$ est une solution de $z' = -2 -2z$ et pour finir on trouve que la solution maximale de l'équation de départ est $y(t)= \\frac{1}{\\sqrt{2e^{-2t}-1}}$ pour $t\\in]-\\infty,\\log(\\sqrt{2})[$.\n",
    "\n",
    "Essayons de représenter ceci graphiquement, comme dans les exercices précédents."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": false,
    "editable": false
   },
   "source": [
    "1) Définir la fonction `f` de l'équation de Cauchy."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "deletable": false
   },
   "outputs": [],
   "source": [
    "# Définition de la fonction f(y,t) de l'équation y'=f(y,t)\n",
    "\n",
    "# test\n",
    "assert f(0, 0) == 0, \"f(0,*) = 0\"\n",
    "assert f(1, -1) == 2, \"f(1,*) = 1\"\n",
    "assert np.shape(f(np.ones(3), np.zeros(3))) == (3, ), \"la fonction doit accepter des vecteurs comme paramètres\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": false,
    "editable": false
   },
   "source": [
    "2) Définir la solution `s` qui vérifie la condition initiale $y(0)=1$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "deletable": false
   },
   "outputs": [],
   "source": [
    "# Définition de la solution s(t) qui vérifie s(0)=1\n",
    "\n",
    "# test\n",
    "assert s(0) == 1, \"la condition initiale est s(0) = 1\"\n",
    "assert np.shape(s(np.zeros(3))) == (3, ), \"la fonction doit accepter des vecteurs comme paramètres\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": false,
    "editable": false
   },
   "source": [
    "3) Adapter les représentations graphiques dynamiques des exercices précédents à cet exercice, représentant la solution `s`, ainsi que la solution trouvée par `odeint` qui vérifie $y(t_0) = y_0$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "deletable": false
   },
   "outputs": [],
   "source": [
    "# choix de la fenêtre\n",
    "tmin, tmax, ymin, ymax = -2, 2, -2, 2\n",
    "\n",
    "# choix de la condition initiale grâce à des curseurs\n",
    "@interact(t_0=(tmin, tmax, .5), y_0=(ymin, ymax, .25))\n",
    "# le dessin pour la condition initiale (t_0,y_0)\n",
    "def dyn_solution(t_0=0, y_0=0):\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": false,
    "editable": false
   },
   "source": [
    "## Amélioration \n",
    "\n",
    "On constate que les solutions trouvées par `odeint` explosent numériquement. Une façon « d'arranger » ça est de « ralentir » `f` en dehors de la fenêtre d'affichage. Pour cela on définit une fonction `ralentir` qui modifie `f` en dehors de la fenêtre visible pour éviter les explosions en temps « longs »."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "deletable": false,
    "editable": false
   },
   "outputs": [],
   "source": [
    "# cutoff une fonction C^0 qui :\n",
    "#  - est égale à 1 pour |Y| < r\n",
    "#  - décroit exponentiellement en dehors de ce disque\n",
    "def cutoff(Y, r):\n",
    "    nY = np.maximum(Y * Y - r * r, 0)\n",
    "    return np.exp(-nY)\n",
    "\n",
    "# la fonction qui « ralentie » F en dehors de l'intervalle [ymin, ymax]\n",
    "def ralentir(f, ymin, ymax):\n",
    "    r = max(abs(ymax), abs(ymin))  # [-r,r] contient [ymin,ymax]\n",
    "    return lambda Y, t: f(Y,t) * cutoff(Y, r)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": false,
    "editable": false
   },
   "source": [
    "4) Reprendre le code précédent, mais en remplaçant `f` par `F = ralentir(f, ymin, ymax)` dans `odeint`. Ainsi on devrait éviter les explosions numériques."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "deletable": false
   },
   "outputs": [],
   "source": [
    "# choix de la fenêtre\n",
    "tmin, tmax, ymin, ymax = -2, 2, -2, 2\n",
    "\n",
    "# choix de la condition initiale grâce à des curseurs\n",
    "@interact(t_0=(tmin, tmax, .5), y_0=(ymin, ymax, .25))\n",
    "# le dessin pour la condition initiale (t_0,y_0)\n",
    "def dyn_solution(t_0=0, y_0=0):\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
}