{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "\\newcommand{\\RR}{\\mathbb{R}}\n",
    "\\newcommand{\\R}{\\mathbb{R}}\n",
    "\\newcommand{\\N}{\\mathbb{N}}\n",
    "$$\n",
    "\n",
    "\n",
    "## L3 Mathématiques - 2023/24 - Univ. Paris-Saclay\n",
    "## Analyse théorique et numérique des EDO \n",
    "---------------------------------------------------------------------------------------------------------------\n",
    "\n",
    "## TP 3 : schémas numériques pour les EDO. Convergence et ordre.\n",
    "---------------------------------------------------------------------------------------------------------------\n",
    "\n",
    "## Partie 1 - CC. \n",
    "#### <font color='red'>A déposer à la fin sur ecampus, dans la page principal de MEU304, dans le dépôt nommé CC2 groupe_mardi_15h45 (ou à envoyer par mail à jean-baptiste.apoung-kamga@universite-paris-saclay.fr ).</font>\n",
    "\n",
    "\n",
    "Soient $t_0,\\ y_0\\in\\R,\\ f:\\R\\times\\R\\longrightarrow\\R$, et considérons le problème de Cauchy  \n",
    "\\begin{equation*}\n",
    "(PC)\\ \\ \\ \\ \\ \\begin{cases}\n",
    "  y'(t) & =  f(t,y(t)), \\\\\n",
    "  y(t_0) & =  y_0.\n",
    "\\end{cases} \n",
    "\\end{equation*}\n",
    "\n",
    "On suppose que ce problème a une unique solution $y:J\\longrightarrow\\mathbb{R}$ définie dans un intervalle $J\\subseteq\\mathbb{R}.$ On s'intéresse à approcher $y$ dans un intervalle de la forme $[t_0, t_0+T]\\subseteq J$, avec $T>0$. \n",
    "\n",
    "Pour ce faire : \n",
    "\n",
    " - on se donne $N\\in\\mathbb{N}$ ;\n",
    " \n",
    " - on définit $h=\\frac{T}{N}$  ;\n",
    " \n",
    " - on définit une sub-division de $[t_0,t_0+T]$ donnée par les **N+1** points $t_n=t_0+n h,\\ n=0,\\dots,N$ ;\n",
    "\n",
    " - on définit des valeurs approchées $y^0,\\ y^1,\\cdots,\\ y^N$ de la solution exacte $t\\mapsto y(t)$ de $(PC)$ aux points $t_0,\\ t_1,\\cdots,\\ t_N$ par le schéma suivant :\n",
    "\n",
    "$$\n",
    "\\begin{cases}\n",
    "y^0&=y_0,\\\\\n",
    "y^{n+1}&=y^{n} + \\frac h2 \\big(f(t_n,y^{n})+f(t_{n+1},\\tilde{y}_{n+1})\\big),\\\\\n",
    "&\\ \\ \\ \\ \\textrm{où}\\ \\ \\ \\tilde{y}_{n+1}=y^{n} + h f\\big(t_n+\\frac h2,y^{n}+\\frac h2 f(t_{n},y^n)\\big),\\ \\ \\ \\ \\ \\ \\textrm{pour }\\ n\\in\\{0,\\dots,N-1\\}.\n",
    "\\end{cases}\n",
    "$$  \n",
    "\n",
    "**Question 1.** Écrire une fonction python de la forme\n",
    "\n",
    "`mon_schema(t0, T, y0, h, f)`\n",
    "\n",
    "(ou de la forme `mon_schema(t0,T,y0,N,f)`, si vous préférez que votre fonction prenne en argument $N$ au lieu de $h$),\n",
    "prenant en argument  l'instant initiale **t0**, la longueur **T** de l'intervalle de temps $[t_0,t_0+T],$ la donnée initiale $y_0$, le pas de temps $h$ (ou le nombre de points $N$) et la fonction $f$ définissant le second membre de l'EDO. Cette fonction devra retourner deux vecteurs : \n",
    "\n",
    "- $[t_0,\\, t_1,\\, ...,\\, t_N],$  tableau `numpy` unidimensionnel de taille $(N+1)$ représentant la subdivision de l'intervalle $[t_0,t_0+T]$ de pas $h$ considérée, \n",
    "\n",
    "- $[y_0,\\, y_1,\\, ...,\\, y_N],$ tableau `numpy` de taille $(N+1)$ représentant la solution approchée aux instants $t_n,\\ n=0,\\dots,N.$\n",
    "\n",
    "(vous pouvez aussi définir ces vecteurs sous forme de listes).\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "## Question 1\n",
    "def mon_schema(t0, T, y0, h, f):\n",
    "    t = np.arange(t0, t0 + T + h, h)\n",
    "    y = 0 * t\n",
    "    y[0] = y0\n",
    "    for n in range(len(t)-1):\n",
    "        yn1 = y[n] + h * f(t[n] + h / 2, y[n] + h/2 * f(t[n], y[n]))\n",
    "        y[n+1] = y[n] + h / 2 * (f(t[n], y[n]) + f(t[n+1], yn1))\n",
    "    return t, y"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Question 2.** On considère le problème de Cauchy \n",
    "\n",
    "\\begin{equation*}\n",
    "(P)\\ \\ \\ \\ \\ \\begin{cases}\n",
    "  y'(t) & =  y+\\sin(t)y^2 \\\\\n",
    "  y(0) & =  1.\n",
    "\\end{cases} \n",
    "\\end{equation*}\n",
    "\n",
    "On souhaite approcher la solution de cette équation **dans l'intervalle $[0,0.5]$**.\n",
    "\n",
    "*Remarque : la solution exacte de $(P)$ est $y(t)=\\displaystyle{\\frac{1}{\\frac12e^{-t}+\\frac{\\cos(t)-\\sin(t)}{2}}}.$*\n",
    "\n",
    "\n",
    "**Question 2.** Calculer la solution approchée de $(P)$ obtenue par cette méthode numérique dans l'intervalle de temps $[0,0.5]$, pour $N=100$. Tracer sur une même fenêtre graphique la solution exacte et la solution approchée. Rajouter des légendes et un titre à votre graphique."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x7f39e37a3b50>"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "[encoded data, 36396 characters]",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "## Question 2\n",
    "\n",
    "# f second membre de l'EDO\n",
    "def f(t,y):\n",
    "    return y + np.sin(t) * np.power(y, 2)\n",
    "\n",
    "# sol. exacte de (P)\n",
    "def yex(t):\n",
    "    return 1 / (1/2 * np.exp(-t) + (np.cos(t)-np.sin(t)) / 2)\n",
    "\n",
    "t0, T = 0, 0.5\n",
    "N = 100\n",
    "y0 = 1\n",
    "t, y_app = mon_schema(t0, T, y0, T/N, f)\n",
    "y_ex = yex(t)\n",
    "\n",
    "plt.figure(1)\n",
    "plt.plot(t, y_ex, label='Solution exacte', lw=3)\n",
    "plt.plot(t, y_app, '+ ', label=f'Solution approchée pour N={N}')\n",
    "plt.title('Solutions approchées pour mon_schema')\n",
    "plt.xlabel(f'$t$')\n",
    "plt.ylabel(f'$y$')\n",
    "plt.legend(fontsize=7)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Question 3.** Calculez les solutions approchées de $(P)$ obtenues avec cette méthode numérique, avec $h=0.5/2^s$ (ou $N=2^{s}$) pour $s = 2,3,...,10$. \n",
    "\n",
    "- **3.1** Représentez dans la même figure la différence en valeur absolue entre la solution exacte et la solution approchée, en fonction du temps, pour chaque valeur de $h$. \n",
    "\n",
    "- **3.2** Calculez, pour chaque valeur de $h$, l'erreur globale \n",
    "\\begin{equation*}\n",
    "E_h =  \\max_{n=0,\\cdots,N}( |y(t_n)-y^{n}|), \n",
    "\\end{equation*} avec $t\\mapsto y(t)$ la solution exacte de $(P)$, et, pour chaque $h=T/N$, $(y^0,\\dots,y^N)$ la solution approchée associée à une discrétisation de pas $h$.  \n",
    "- **3.2.1** Représentez dans la même figure $E_h/h$ et $E_h/h^2$ en fonction de $N$ (pour ce faire, gardez les valeurs de $E_h$, de $h$ et de $N$, pour $s=2,...,10$, dans des listes ou des tableaux numpy).\n",
    "- **3.2.2** Représentez dans  une autre figure, en échelle logarithmique, l'erreur $E_h$ en fonction du pas de temps $h$, autrement dit, représentez $\\log(E_h)$ en fonction de $\\log(h)$. Vous devez obtenir des points qui sont à peu près alignés sur une droites. Vérifiez graphiquement que c'est le cas, en représentant des droites de pente $1$ et $2$ (pour ce faire, représentez respectivement $log(h)$ en fonction de $log(h)$ et $2log(h)$ en fonction de $log(h)$). Vérifiez que tous les points semblent être dans une droite paralèlle à une de ces deux droites.\n",
    "\n",
    "- **3.3** Que semblent montrer les figures des questions 3.2.1 et 3.2.2 sur la convergence du schéma ?\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0.5, 1.0, 'Solutions approchées de (P) obtenus avec mon_schema pour différentes valeurs du pas h')"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "[encoded data, 59184 characters]",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Question 3\n",
    "plt.figure(2)\n",
    "\n",
    "for s in range(2,11):\n",
    "    h=(1/2)/(2**s)\n",
    "    t, y_app = mon_schema(t0, T, y0, h, f)\n",
    "    plt.plot(t, y_app, label=f'h={h}')\n",
    "\n",
    "plt.legend(fontsize=7)\n",
    "plt.xlabel(f'$t$')\n",
    "plt.ylabel('$y^n$')\n",
    "plt.title('Solutions approchées de (P) obtenus avec mon_schema pour différentes valeurs du pas h')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0.5, 1.0, 'différence en valeur absolue entre sol. exacte et sol. approchée par mon_schema, pour différentes valeurs du pas h')"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "[encoded data, 74936 characters]",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "E=[]\n",
    "H=[]\n",
    "\n",
    "plt.figure(3)\n",
    "for s in range(2,11):\n",
    "    h=(1/2)/(2**s)\n",
    "    t, y_app = mon_schema(t0, T, y0, h, f)\n",
    "    y_ex = yex(t)\n",
    "    err = np.max(np.abs(y_app - y_ex))\n",
    "    plt.plot(t, np.abs(y_app - y_ex), label=f'h={h}')\n",
    "    E.append(err)\n",
    "    H.append(h)\n",
    "    \n",
    "plt.legend()\n",
    "plt.xlabel(f'$t$')\n",
    "plt.ylabel('$|y(t_n) - y^n|$')\n",
    "plt.title('différence en valeur absolue entre sol. exacte et sol. approchée par mon_schema, pour différentes valeurs du pas h')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0.5, 1.0, 'Erreur globale en fonction de N')"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "[encoded data, 26724 characters]",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure(4)\n",
    "\n",
    "H, E = np.array(H), np.array(E)\n",
    "\n",
    "plt.plot(T/H, E/H, label=f'$E_h / h$')\n",
    "plt.plot(T/H, E/H**2, label=f'$E_h / h^2$')\n",
    "\n",
    "plt.legend()\n",
    "plt.xlabel('N')\n",
    "plt.ylabel('Erreur globale')\n",
    "plt.title('Erreur globale en fonction de N')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0, 0.5, '$log(E)$')"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "[encoded data, 58788 characters]",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure(5)\n",
    "\n",
    "plt.plot(np.log(H), np.log(E), '+-', label='erreur')\n",
    "plt.plot(np.log(H), np.log(H), '-', label='droite pente 1')\n",
    "plt.plot(np.log(H), 2*np.log(H), '-', label='droite pente 2')\n",
    "\n",
    "plt.legend()\n",
    "plt.title('Erreur pour la méthode mon_schema en echelle logarithmique : log(E) en fonction de log(h)')\n",
    "plt.xlabel(f'$log(h)$')\n",
    "plt.ylabel(f'$log(E)$')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Question 3.3\n",
    "\n",
    "Avec les schemas 3.2.1 et 3.2.2, on remarque que notre schema est d'ordre 2.\n",
    "En effet, dans le schema 3.2.1, l'erreur globale d'ordre 1 est nulle quand N augmente alors que l'erreur globale d'ordre 2 est environ égale a 1.5 quand N augmente.\n",
    "De plus, dans le schema 3.2.2, la droite d'erreur est parallèle à la droite de pente 2 donc notre schema est convergent d'ordre 2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Partie 2. Comparaison de différentes méthodes numériques.\n",
    "\n",
    "L'objectif de cette partie est de mettre en oeuvre différents schémas numériques pour approcher la solution d'une équation différentielle ordinaire, et de comprendre des notions comme l'ordre de précision d'une méthode numérique. \n",
    "\n",
    "**Quelques conseils pour l'organisation de vos programmes :** *vous allez programmer plusieurs schémas numériques et les tester sur plusieurs exemples d'EDO. Il est souhaitable de bien organiser votre programme pour ne pas se perdre. Par exemple vous pourrez définir au tout début les fonctions second membre des EDO que vous allez tester, en les appelant par exemple $f1,\\ f2,\\dots.$ Vos fonctions doivent avoir toujours deux variables, $t$ et $y$, même si elles ne dépendent pas de $t$ (comme ça les fonctions définissant les schémas seront adaptées à tous les cas).Vous pourrez ensuite définir les solutions exactes de ces équations, si elles sont connues, en les appelant par exemple $yex1,\\ yex2,\\dots$. Également, vous pouvez définir ensuite les fonctions que vous allez construire pour chacun des schémas. Lorsque vous faites appelle à une fonction pour tester un certain schéma dans une situation concrète, il faut définir avant les paramètres du cas que vous allez tester. Par exemple, faire*\n",
    "\n",
    "`t0=0`\n",
    "\n",
    "`tf=1`\n",
    "\n",
    "`y0=2`\n",
    "\n",
    "`h=0.5`\n",
    "\n",
    "`euler_exp(f1, t0, tf, y0, h)`\n",
    "\n",
    "*au lieu de*\n",
    "\n",
    "`euler_exp(f1,0,1,2,0.5)`.\n",
    "\n",
    "*Cela vous permet de facilement changer les paramètres s'il le faut, et la visibilité de votre programme sera aussi facilitée.*\n",
    "\n",
    "On se donne alors le problème de Cauchy  \n",
    "\\begin{equation*}\n",
    "(PC)\\ \\ \\ \\ \\ \\begin{cases}\n",
    "  y'(t) & =  f(t,y(t)), \\\\\n",
    "  y(t_0) & =  y_0\n",
    "\\end{cases} \n",
    "\\end{equation*}\n",
    "\n",
    "où $f:\\R\\times\\R^n\\longrightarrow\\R^n$ (si $n>1$, il s'agit d'un système d'équations différentielles).\n",
    "\n",
    "Avec les mêmes notations de la partie 1, on va considérer les schémas numériques suivants pour approcher la solution exacte de ce problème :\n",
    "\n",
    "-  **Euler explicite** : $y^{n+1}=y^{n}+h f(t_n,y^{n})$ ; \n",
    "  \n",
    "-  **Euler implicite** : $y^{n+1}=y^{n}+h f(t_{n+1},y^{n+1})$  ;\n",
    "\n",
    "-  **Heun** :  $y^{n+1}=y^{n}+\\frac{h}{2}\\big(p_1+p_2\\big)$, où\n",
    "  $p_1=f(t_n,y^{n})$, $\\ p_2=f(t_n+h,y^{n}+h p_1)$ ;\n",
    "\n",
    "- **Crank-Nicolson** : $y^{n+1}=y^{n}+\\frac{h}{2}\\big(f(t_{n},y^{n})+f(t_{n+1},y^{n+1})\\big)$ ;\n",
    "- **Runge Kutta 4** : $y^{n+1}=y^{n}+\\frac{h}{6}(p_1+2p_2+2p_3+p_4)$, où $p_1=f(t_n,y^{n})$, $\\ p_2=f(t_n+\\frac{h}{2},y^{n}+\\frac{h}{2} p_1)$, $\\ p_3=f(t_n+\\frac{h}{2},y^{n}+\\frac{h}{2} p_2)$ et \n",
    "$\\ p_4=f(t_n+h,y^{n}+h p_3)$. \n",
    "- **Adams-Bashforth à deux pas** : $y^{n+1}=y^{n}+\\frac{h}{2}\\big(3f(t_n,y^{n})-f(t_{n-1},y^{n-1})\\big)$.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercice 1.** \n",
    "\n",
    "Écrire, pour chaque  schéma explicite à un pas, une fonction python de la forme\n",
    "\n",
    "`nom_du_schema(t0, T, y0, h, fct)`\n",
    "\n",
    "(ou de la forme `mon_schema(t0,T,y0,N,fct)`),\n",
    "prenant en argument  l'instant initiale **t0**, la longueur **T** de l'intervalle de temps $[t_0,t_0+T],$ la donnée initiale $y_0$, le pas de temps $h$ (ou le nombre de points $N$) et la fonction `fct` définissant le second membre de l'EDO. Cette fonction devra retourner deux tableaux :  \n",
    "- $[t_0,\\, t_1,\\, ...,\\, t_N],$ tableau `numpy` unidimensionnel de taille $N+1$ représentant la subdivision de l'intervalle $[t_0,t_0+T]$ de pas $h$ considérée, \n",
    "- $[y^0,\\, y^1,\\, ...,\\, y^N],$ tableau `numpy` de taille $(N+1)\\times n$ représentant la solution approchée aux instants $t_n,\\ n=0,\\dots,N.$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "def euler_implicite(t0, T, y0, h, f):\n",
    "    t = np.arange(t0, t0 + T + h, h)\n",
    "    y = 0 * t\n",
    "    y[0] = y0\n",
    "    for n in range(len(t)-1):\n",
    "        y[n+1] = y[n] + h * f(t[n+1], y[n+1])\n",
    "    return t, y\n",
    "\n",
    "def crank(t0, T, y0, h, f):\n",
    "    t = np.arange(t0, t0 + T + h, h)\n",
    "    y = 0 * t\n",
    "    y[0] = y0\n",
    "    for n in range(len(t)-1):\n",
    "        y[n+1] = y[n] + h/2 * (f(t[n], y[n]) + f(t[n+1], y[n+1]))\n",
    "    return t, y\n",
    "\n",
    "def adams(t0, T, y0, h, f):\n",
    "    t = np.arange(t0, t0 + T + h, h)\n",
    "    y = 0 * t\n",
    "    y[0] = y0\n",
    "    for n in range(1, len(t)-1):\n",
    "        y[n+1] = y[n] + h/2 * (3* f(t[n], y[n]) - f(t[n-1], y[n-1]))\n",
    "    return t, y"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "def euler_explicite(t0, T, y0, h, f):\n",
    "    N = int(T/h)\n",
    "    n = len(y0)\n",
    "    t = np.linspace(t0, t0 + T, N+1)\n",
    "    y = np.zeros((N+1, n))\n",
    "    y[0,] = y0\n",
    "    for n in range(N):\n",
    "        y[n+1] = y[n] + h * f(t[n], y[n])\n",
    "    return t, y\n",
    "\n",
    "def heun(t0, T, y0, h, f):\n",
    "    N = int(T/h)\n",
    "    n = len(y0)\n",
    "    t = np.linspace(t0, t0 + T, N+1)\n",
    "    y = np.zeros((N+1, n))\n",
    "    y[0,] = y0\n",
    "    for n in range(N):\n",
    "        p1 = f(t[n], y[n])\n",
    "        p2 = f(t[n] + h, y[n] + h * p1)\n",
    "        y[n+1] = y[n] + h/2 * (p1 + p2)\n",
    "    return t, y\n",
    "\n",
    "def runge(t0, T, y0, h, f):\n",
    "    N = int(T/h)\n",
    "    n = len(y0)\n",
    "    t = np.linspace(t0, t0 + T, N+1)\n",
    "    y = np.zeros((N+1, n))\n",
    "    y[0,] = y0\n",
    "    for n in range(N):\n",
    "        p1 = f(t[n], y[n])\n",
    "        p2 = f(t[n] + h/2, y[n] + h/2 * p1)\n",
    "        p3 = f(t[n] + h/2, y[n] + h/2 * p2)\n",
    "        p4 = f(t[n] + h, y[n] + h* p3)\n",
    "        y[n+1] = y[n] + h/6 * (p1 + 2*p2 + 2*p3 + p4)\n",
    "    return t, y"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    " **Exercice 2.** \n",
    "\n",
    "**Question 1.** Testez les schémas d'Euler explicite, d'Heun et RK4 sur le modèle logistique\n",
    "$$\n",
    "(P_1)\\ \\ \\ \\ \n",
    "\\begin{cases}\n",
    "y'(t)=c y (1 - \\frac{y}{b}),\\\\\n",
    "y(0)=a,\n",
    "\\end{cases}\n",
    "$$\n",
    "dont la solution exacte est\n",
    "$$\n",
    "y(t) = \\frac{b}{1 + \\frac{b-a}{a} e^{-ct}},\n",
    "$$\n",
    "avec les données $c=1$, $b=2$, $a=0.1$, dans l'intervalle $[0,15]$, avec un pas $h=0.2$. Tracer sur la même fenêtre la solution exacte et la solution approchée obtenue avec chacun des schémas. Tracer dans une autre figure, pour chacun des schémas, la différence en valeur absolue entre la solution exacte et la solution approchée. Quelle méthode semble être la meilleure ?\n",
    "\n",
    "\n",
    "**Question 2.** Testez ensuite les 3 schémas dans le cas vectoriel $n>1$ ($n=2$) sur le problème\n",
    "$$\n",
    "(P_2)\\ \\ \\ \\\n",
    "\\begin{cases}\n",
    "y''(t)  =    -y(t) + \\cos(t)  \\\\\n",
    "y(0)  =   5,\\ \\ y'(0)  =   1,  \n",
    "\\end{cases}\n",
    "$$\n",
    "dont la solution exacte est \n",
    "$$\n",
    "y(t) = \\frac{1}{2} \\sin(t) t  + 5 \\cos(t) + \\sin(t),\n",
    "$$\n",
    "dans l'intervalle $[0,15]$, avec un pas  $h=0.2$.\n",
    "\n",
    "Pour ce faire, il faut écrire l'équation d'ordre 2 de $(P_2)$ sous la forme d'un système de deux équations d'ordre 1 dans les nouvelles variables $u(t)=y(t)$ et $v(t)=y'(t)$. On se ramènera alors à la résolution d'une équation de la forme\n",
    "$$\n",
    "X'=F(t,X),\\ \\ \\ \\textrm{avec}\\ X=(u,v)=(y,y')^T.\n",
    "$$\n",
    "Représenter à nouveau la solution exacte et la solution approchée dans une même fenêtre graphique.\n",
    "\n",
    "*La solution $y$ de $(P_2)$ correspond à la première composante du vecteur $X$ ci-dessus. Votre fonction `schema` retournera dans ce cas, si vous avez respecté la structure conseillée, un tableau de taille $(N+1)\\times2$, $N$ étant le nombre de points de la discrétisation. Ce tableau donne les valeurs approchées de $X$ au points de la discrétisation considérée, la première colonne correspondant à la première composante de $X$, la seconde à la seconde composante de $X$. La solution approchée de $(P_2)$ que l'on cherche correspond alors à la première colonne de ce tableau.*\n",
    "\n",
    "**Question 3.** Étude de l'erreur et convergence.\n",
    "\n",
    "Considérons le problème\n",
    "$$\n",
    "(P_3)\\ \\ \\ \\\n",
    "\\begin{cases}\n",
    "y'(t)=\\frac{\\cos(t)-y(t)}{1+t},\\\\\n",
    "y(0)=-\\frac14, \n",
    "\\end{cases}\n",
    "$$\n",
    "dont la solution exacte est \n",
    "$$\n",
    "y(t) = \\frac{\\sin(t)-1/4}{1 + t}.\n",
    "$$\n",
    "\n",
    "Pour chacune des 3 méthodes (EE, Heun et RK4) :\n",
    "\n",
    "**3.1** Calculez les solutions approchées de $(P_3)$ obtenues avec avec $h=1/2^s$ pour $s = 2,...,10$, sur l'intervalle de temps $[0,10]$. \n",
    "\n",
    "**3.2** Calculer, pour chaque valeur de $h=1/2^s,\\ s = 1,2,...,10,$ l'erreur globale $E_h$.  \n",
    "Représentez ensuite dans une même figure, l'erreur en fonction du pas de temps $h$ en échelle logarithmique  (autrement dit, représentez $\\log(E_h)$ en fonction de $\\log(h)$), pour chaque méthode. Vous devez obtenir des points qui sont à peu près alignés sur des droites. Vérifiez graphiquement que c'est le cas, en estimant la pente de la droite passant au plus prêt des points (ou en représentant des droites de pente $1,\\ 2$ et $4$ qui passent par un des points et en vérifiant que tous les points sont à peu près sur ces droites).\n",
    "\n",
    "*Si les points $(log(h),log(E_h))$ sont alignés sur une droite de pente $p$, ceci signifie que $\\log(E_h) \\sim C+p\\log(h)=C+\\log(h^p)$ et donc que $E_h\\sim \\widetilde{C}h^p$, pour certaines constantes $C$ et $\\widetilde{C}$. On dit que la méthode est d'ordre p : c'est l'ordre de la puissance de $h^p$ dans cette relation. On a donc que l'erreur globale $E_h$ tend vers 0 comme $h^p$. L'ordre d'une méthode donne une indication sur sa vitesse de convergence. Une méthode d'ordre $p$ est une méthode dont l'erreur globale tend vers $0$ comme $h^p$. Donc plus l'ordre est élevé, plus la méthode converge plus vite*.\n",
    "\n",
    "\n",
    "**Remarque** : *pour étudier numériquement l'ordre d'une méthode, on utilise souvent l'échelle logarithmique pour tracer l'erreur en fonction du pas de discrétisation $h$. La pente de la droite obtenue donne l'ordre $p$ de la méthode : si $E_h \\sim Ch^p$ alors $\\log(E_h)\\sim \\log(C) + p\\log(h)$.*\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0.5, 1.0, 'différence en valeur absolue entre sol. exacte et sol. approchée, pour différents schemas')"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "[encoded data, 48988 characters]",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 50268 characters]",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Question 1\n",
    "a, b, c = 0.1, 2, 1\n",
    "\n",
    "# f second membre de l'EDO\n",
    "def f1(t,y):\n",
    "    return c * y * (1 - y/b)\n",
    "\n",
    "# sol. exacte de (P1)\n",
    "def yex1(t):\n",
    "    return b / (1 + (b-a)/a * np.exp(-c*t))\n",
    "\n",
    "t0, T = 0, 15\n",
    "h = 0.2\n",
    "y0 = np.array([a])\n",
    "\n",
    "plt.figure()\n",
    "for schema in [euler_explicite, heun, runge]:\n",
    "    t, y_app = schema(t0, T, y0, h, f1)\n",
    "    plt.plot(t, y_app.ravel(), '+ ', label=f'Solution approchée par {schema.__name__}')\n",
    "\n",
    "t = np.arange(t0, t0+T+h, h)\n",
    "y_ex = yex1(t)\n",
    "plt.plot(t, y_ex, label='Solution exacte', lw=1)\n",
    "plt.title(f'Solutions approchées pour le schema {schema.__name__}')\n",
    "plt.xlabel(f'$t$')\n",
    "plt.ylabel(f'$y$')\n",
    "plt.legend(fontsize=7)\n",
    "\n",
    "plt.figure()\n",
    "for schema in [euler_explicite, heun, runge]:    \n",
    "    t, y_app = schema(t0, T, y0, h, f1)\n",
    "    plt.plot(t, np.abs(y_app.ravel() - yex1(t)), label=f'Schema {schema.__name__}')\n",
    "    \n",
    "plt.legend()\n",
    "plt.xlabel(f'$t$')\n",
    "plt.ylabel('$|y(t_n) - y^n|$')\n",
    "plt.title(f'différence en valeur absolue entre sol. exacte et sol. approchée, pour différents schemas')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x7f39d95d84f0>"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "[encoded data, 52596 characters]",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Question 2\n",
    "def F(t, Y):\n",
    "    Z = np.zeros_like(Y)\n",
    "    Z[0] = Y[1]\n",
    "    Z[1] = -Y[0] + np.cos(t)\n",
    "    return Z\n",
    "\n",
    "def yexF(t):\n",
    "    return 1/2 * np.sin(t) * t + 5 * np.cos(t) + np.sin(t),\n",
    "\n",
    "t0, T = 0, 15\n",
    "h = 0.2\n",
    "Y0 = np.array([5,1])\n",
    "\n",
    "plt.figure()\n",
    "for schema in [euler_explicite, heun, runge]:\n",
    "    t, y_app = schema(t0, T, Y0, h, F)\n",
    "    plt.plot(t, y_app[:,0], '+ ', label=f'Solution approchée par {schema.__name__}')\n",
    "\n",
    "t = np.arange(t0, t0+T+h, h)\n",
    "y_ex = yexF(t)\n",
    "plt.plot(t, y_ex[0], label='Solution exacte', lw=3)\n",
    "plt.title('Solutions approchées par differents schemas')\n",
    "plt.xlabel(f'$t$')\n",
    "plt.ylabel(f'$y$')\n",
    "plt.legend(fontsize=7)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "# Question 3\n",
    "\n",
    "# f second membre de l'EDO\n",
    "def f3(t,y):\n",
    "    return (np.cos(t) - y) / (1+t)\n",
    "\n",
    "# sol. exacte de (P1)\n",
    "def yex3(t):\n",
    "    return (np.sin(t) - 1/4) / (1 + t)\n",
    "\n",
    "t0, T = 0, 10\n",
    "y0 = np.array([-1/4])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<Figure size 640x480 with 0 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 88992 characters]",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 59204 characters]",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 59340 characters]",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure()\n",
    "for schema in [euler_explicite, heun, runge]:\n",
    "    plt.figure()\n",
    "    for s in range(2,11):\n",
    "        h=1/(2**s)\n",
    "        t, y_app = schema(t0, T, y0, h, f3)\n",
    "        plt.plot(t, y_app, label=f'h={h}')\n",
    "\n",
    "    plt.legend(fontsize=7)\n",
    "    plt.xlabel(f'$t$')\n",
    "    plt.ylabel('$y^n$')\n",
    "    plt.title(f'Solutions approchées de (P) obtenus avec {schema.__name__} pour différentes valeurs du pas h')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 66364 characters]",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 67132 characters]",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 57532 characters]",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "for schema in [euler_explicite, heun, runge]:\n",
    "    H, E = [], []\n",
    "    plt.figure()\n",
    "    for s in range(2,11):\n",
    "        h=1/(2**s)\n",
    "        t, y_app = schema(t0, T, y0, h, f3)\n",
    "        E.append(np.max(np.abs(yex3(t) - y_app.ravel())))\n",
    "        H.append(h)\n",
    "    \n",
    "    H, E = np.array(H), np.array(E)\n",
    "    plt.plot(np.log(H), np.log(E), '+-', label='erreur')\n",
    "    plt.plot(np.log(H), np.log(H), '-', label='droite pente 1')\n",
    "    plt.plot(np.log(H), 2*np.log(H), '-', label='droite pente 2')\n",
    "    plt.plot(np.log(H), 3*np.log(H), '-', label='droite pente 3')\n",
    "    \n",
    "    plt.legend()\n",
    "    plt.title(f'Erreur pour la méthode {schema.__name__} en echelle logarithmique : log(E) en fonction de log(h)')\n",
    "    plt.xlabel(f'$log(h)$')\n",
    "    plt.ylabel(f'$log(E)$')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercice 3. Problèmes raides**\n",
    "\n",
    "Nous nous intéressons ici au cas du modèle de Van der Pol :\n",
    "\n",
    "- Le modèle de Van Der Pol \n",
    "$$\n",
    "(P)\\ \\ \\ \\ \n",
    "\\begin{cases}\n",
    "y''(t)  =    (a-y(t)^2) y'(t) - y(t),\\\\\n",
    "y(0)  =   \\frac{1}{10},\\ \\ y'(0)  = 1,  \n",
    "\\end{cases}\n",
    "$$\n",
    "dont la solution analytique on ne sait pas calculer, même si elle existe (c'est une conséquence du théorème de Cauchy-Lipschitz). On prendra $a=5$ ou $a=8$, $[t_0,t_f] =[0, 20]$ dans les tests numériques.\n",
    "\n",
    "**Remarque :** *il n'est pas simple de définir ce que c'est une équation raide. Il s'agit d'équations dont les solutions présentent de grandes variations initiales ou de systèmes d'équations où les échelles de variation des différentes inconnues sont très différentes. Il s'agit d'équations qui sont en général difficiles d'approcher avec des schémas explicites d'ordre bas.*\n",
    "\n",
    "**Question 1.** Écrivez $(P)$ sous la forme d'un système différentiel de deux équations  d'ordre 1. \n",
    "\n",
    "**Question 2.** Testez tous les schémas explicites proposés lorsque l'on prend $a=8$, $t \\in [0, 20]$ et $h=0.1,\\, 0.05,\\, 0.025$ ? Recommencer avec $a=5$. Commentez les résultats obtenus.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/tmp/ipykernel_536/3007116115.py:7: RuntimeWarning: overflow encountered in scalar multiply\n",
      "  return np.array([Y[1], (a - Y[0]**2) * Y[1] - Y[0]])\n",
      "/tmp/ipykernel_536/3007116115.py:7: RuntimeWarning: overflow encountered in scalar power\n",
      "  return np.array([Y[1], (a - Y[0]**2) * Y[1] - Y[0]])\n",
      "/tmp/ipykernel_536/3996006130.py:8: RuntimeWarning: invalid value encountered in add\n",
      "  y[n+1] = y[n] + h * f(t[n], y[n])\n",
      "/tmp/ipykernel_536/3996006130.py:20: RuntimeWarning: invalid value encountered in add\n",
      "  y[n+1] = y[n] + h/2 * (p1 + p2)\n"
     ]
    },
    {
     "data": {
      "image/png": "[encoded data, 38616 characters]",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 36648 characters]",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 57456 characters]",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "a = 5\n",
    "t0, T = 0, 20\n",
    "H = [0.1, 0.05, 0.025]\n",
    "Y0 = np.array([1/10, 1])\n",
    "\n",
    "def P(t, Y):\n",
    "    return np.array([Y[1], (a - Y[0]**2) * Y[1] - Y[0]])\n",
    "\n",
    "for schema in [euler_explicite, heun, runge]:\n",
    "    plt.figure()\n",
    "    for h in H:\n",
    "        t, y_app = schema(t0, T, Y0, h, P)\n",
    "        plt.plot(t, y_app[:,0], '+ ', label=f'Solution approchée pour h={h}')\n",
    "    plt.legend()\n",
    "    plt.title(f'Solutions approchees pour differents pas h par {schema.__name__}')\n",
    "    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercice 4. Méthodes multi-pas.** Les méthodes multi-pas, comme la méthode d'Adams-Bashforth, utilisent d'autres points que $(t^n,y(t^n))$ pour calculer une valeur approchée de $y(t^{n+1}).$ L'intérêt de ces méthodes est d'obtenir des méthodes d'ordre élévé pour une compléxité inférieure que les méthodes de Runge-Kutta.\n",
    "\n",
    "Reprendre les questions de l'exercice précédent pour le schéma d'Adams-Bashforth. On remarque qu'il faut calculer $y^1$ par une autre méthode, on utilisera alors la méthode de Heun. Essayer d'établir numériquement l'ordre de la méthode. \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.10.8"
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {},
   "toc_section_display": true,
   "toc_window_display": false
  },
  "varInspector": {
   "cols": {
    "lenName": 16,
    "lenType": 16,
    "lenVar": 40
   },
   "kernels_config": {
    "python": {
     "delete_cmd_postfix": "",
     "delete_cmd_prefix": "del ",
     "library": "var_list.py",
     "varRefreshCmd": "print(var_dic_list())"
    },
    "r": {
     "delete_cmd_postfix": ") ",
     "delete_cmd_prefix": "rm(",
     "library": "var_list.r",
     "varRefreshCmd": "cat(var_dic_list()) "
    }
   },
   "types_to_exclude": [
    "module",
    "function",
    "builtin_function_or_method",
    "instance",
    "_Feature"
   ],
   "window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
