{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "%config InlineBackend.figure_format = 'retina'\n",
    "\n",
    "import numpy as np                       # pour les numpy array\n",
    "import matplotlib.pyplot as plt          # librairie graphique\n",
    "from scipy.integrate import odeint       # seulement odeint"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1. Un modèle logistique avec la capacité d'accueil du milieu dépendant du temps\n",
    "\n",
    "On rappelle que le modèle logistique s'écrit\n",
    "\\begin{equation*} N'(t) = r \\; N(t) \\Bigl( 1 - \\frac{N(t)}{K} \\Bigr),\n",
    "\\label{logistique} \\tag{1}\n",
    "\\end{equation*}\n",
    "où $r$ et $K$ sont deux constantes strictement positives : $r$ est le taux de croissance (naissance - mort) et $K$ est la capacité du milieu. La solution du problème de Cauchy pour l'équation logistique associé à la donnée initiale $ N( 0 ) = N_0 $ est \n",
    "$$\n",
    "N(t) = \\frac{K}{1 + \\bigl(\\frac{K}{N_0}-1\\bigr) \\text{e}^{-rt}}, \\qquad t\\geq0. \\tag{2}\n",
    "$$\n",
    "\n",
    "Il est peu réaliste de supposer que dans un modèle décrivant la dynamique d’une population humaine la capacité d’accueil du milieu reste toujours constante. On va considérer un modèle où la capacité de milieu est affectée par des changements technologiques.\n",
    "On va supposer que de nouvelles technologies augmentent l’efficacité de consommation de ressources et que elles ne se propagent pas instantanément, mais sont plutôt adoptées à un rythme variable.\n",
    "Ces processus d'adoption sont souvent bien modélisés par une courbe logistique. Premièrement, le taux de l'adoption est lente, car une nouvelle technologie a des difficultés pour remplacer une technologie existante. Ensuite, le taux d'adoption augmente, généralement de façon exponentielle jusqu'à ce que les limites physiques ralentissent l'adoption. \n",
    "\n",
    "Nous allons modéliser la taille de la population avec une fonction $ N $ qui satisfait le problème de Cauchy:\n",
    "$$ N'(t) = r_N \\; N(t) \\Bigl( 1 - \\frac{N(t)}{C(t)} \\Bigr), \\quad N( 0 ) = N_0 \\tag{3}$$\n",
    "avec $ C $, la capacité d'accueil du milieu, solution du problème suivant:\n",
    "$$C'(t) =r ( C - K_* ) \\Bigl( 1 - \\frac{ C - K_* }{ K } \\Bigr), \\quad C( t_{fl} ) = K_0 + K_*, \\tag{4} $$\n",
    "où $ r, K_*, K, t_{fl}, K_0 $ sont des paramètres réels.\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**EXERCICE 1**\n",
    "> 1. Montrer que si $ \\tilde K $ est solution de l'équation logistique (1)  vérifiant $ \\tilde K( t_{fl} ) = K_0 $, alors $ C = \\tilde K + K_* $ est solution de (4).\n",
    "> 2. Montrer que si $ P $ est solution de (1) vérifiant $ P( 0 ) = K_0 $, alors la fonction $ \\tilde K $ définie par $ \\tilde K( t ) = P( t - t_{fl} ) $ est solution de (1) vérifiant $ \\tilde K( t_{fl} ) = K_0 $.\n",
    "> 3. En utilisant l'expression (2) de la solution exacte de l'équation (1) et le résultat des questions 1-2, montrer que la solution exacte de (4) est donnée par la formule\n",
    "$$\n",
    "C( t ) = K_* + \\frac{K} {( 1 + ( K / K_0 - 1 ) * e^{ - r ( t - t_{fl} ) )}}.\n",
    "$$\n",
    "\n",
    "On prend les valeurs suivantes des paramètres: $K= 50$, $K_*=50$, $r = 0.1$, $ t_{fl} = 30 $, $ K_0 = 1 $. \n",
    "> 4. Tracer la courbe de $ C $ sur l'intervalle de temps $ [0, T] $ avec $ T = 130 $. Ajouter les lignes horizontales $ C = K_* $, $ C = K + K_* $ et un marqueur au point $ (t_{fl}, K_0+K_*) $, ajouter une légende. En utilisant la figure proposer une interprétation des paramètres $ K $, $ K_* $, $ r $, $ t_{fl} $, $ K_0 $.\n",
    "\n",
    "Pour la suite on prend les valeurs des paramètres $T=130$, $ N_0 = 10 $, $r_N = 0.2$.\n",
    "> 5. Construire la solution approchée du problème (3) à l'aide de la fonction `odeint` et la tracer (avec des marqueurs). \n",
    "> 6. Superposer la courbe de la fonction $t\\mapsto C(t)$. Que remarquez-vous ?\n",
    "> 7. Faire varier la valeur du paramètre $ N_0 $ et commentez les résultats obtenus."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 90752 characters]\n",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {
      "image/png": {
       "height": 480,
       "width": 571
      }
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 100152 characters]\n",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {
      "image/png": {
       "height": 480,
       "width": 571
      }
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Initialisation des variables\n",
    "T = 130\n",
    "t = np.arange(0, T+1)\n",
    "\n",
    "K = 50       # capacité d'accueil maximale du milieu\n",
    "K_star = 50  # capacité minimale pour maintenir l'espèce\n",
    "r = 0.1      # taux de  corissance de la capacité d'accueil du milieu.\n",
    "t_fl = 30    # la find ed la période de formation.\n",
    "K0 = 1       # Valeur initiale de la capacité d'accueil du milieu.\n",
    "\n",
    "def C(t):\n",
    "    \"\"\"\n",
    "    Fonction retournant la solution exacte du problème au temps t\n",
    "    \"\"\"\n",
    "    return K_star + K / (1 + (K / K0 - 1) * np.exp(-r * (t - t_fl)))\n",
    "\n",
    "# On trace le graphique de la solution exacte\n",
    "plt.plot(t, C(t), label=\"C(t)\")\n",
    "plt.hlines(K_star, 0, T, linestyle='dotted', label=\"C = K*\", color='red')\n",
    "plt.hlines(K + K_star, 0, T, linestyle='dotted', label=\"C = K + K*\", color='green')\n",
    "plt.plot(t_fl, K0 + K_star, 'o', label=\"(t_fl, K0 + K*)\")\n",
    "plt.legend()\n",
    "plt.xlim(0, 130)\n",
    "plt.suptitle(\"Courbe de la solution exacte du problème\")\n",
    "plt.xlabel(\"Temps t\")\n",
    "plt.ylabel(\"Population\")\n",
    "plt.show()\n",
    "\n",
    "# Définition des nouveaux paramètres\n",
    "N0 = 10\n",
    "r_N = 0.2\n",
    "\n",
    "def dN(N, t, C_sol):\n",
    "    \"\"\"\n",
    "    Fonction calculant la dérivée de la solution approchée du problème à l'instant t dépendant de N(t) et de C(t)\n",
    "    \"\"\"\n",
    "    return r_N * N * (1 - N / C_sol(t))\n",
    "\n",
    "t = np.linspace(0, T, 200)\n",
    "N_sol = odeint(dN, N0, t, args=(C,)) # On calcule la solution a l'aide de odeint\n",
    "\n",
    "# On trace le graphique de la solution approchée en comparaison à la solution exacte\n",
    "plt.plot(t, N_sol, label=\"Solution approchée\")\n",
    "plt.plot(t, C(t), label=\"Solution exacte\")\n",
    "plt.legend()\n",
    "plt.xlim(0, 130)\n",
    "plt.suptitle(\"Courbe de la solution approchée du problème\")\n",
    "plt.xlabel(\"Temps t\")\n",
    "plt.ylabel(\"Population\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2. L'impact de la pêche dans le modèle Lotka-Volterra\n",
    "\n",
    "\n",
    "L'objectif de cette partie est d'étudier numériquement l'impact de la peche dans le modèle de Lotka-Volterra et d'illustrer le paradoxe observé par d'Anconna (voir la fiche de cours n°1).\n",
    "\n",
    "On se donne donc deux populations $H(t)$ de sardines (les proies) et $P(t)$ de requins (les prédateurs) au temps $t,$ dont l'évolution au cours du temps suit la loi:\n",
    "\n",
    "\\begin{equation*}\n",
    "     \\begin{cases}\n",
    "        H'(t)&=H(t)( a-p-bP(t)),\\\\\n",
    "        P'(t)&=P(t)( -c -p +dH(t)).\n",
    "     \\end{cases}\n",
    "\\end{equation*}\n",
    "\n",
    "où $a,b,c,d,p >0$. Le coefficient $a$ représente le taux de reproduction des sardines en l'abscence de requins, $c$ le taux de mortalité en l'abscence de sardines, $b$ le taux de mortalités des sardines mangés par les requins et $d$ le taux de reproduction des requins en fonction de la quantité de sardines mangées. Nous pêchons une même proportion $p$ de requins et de sardines ce qui se traduit par l'ajout de deux termes négatifs $-pH(t)$ et $-pP(t)$ dans les équations du modèle Lotka-Volterra.\n",
    "\n",
    "On suppose connues les populations de sardines et de requins à l'instant initial $t_0=0$ :\n",
    "\n",
    "\\begin{equation*}\n",
    "     \\begin{cases}\n",
    "    H(0) = H_0, \\\\\n",
    "    P(0) = P_0,\n",
    "     \\end{cases}\n",
    "\\end{equation*}\n",
    "où $H_0\\ge 0$ et $P_0\\ge 0$ sont donnés.\n",
    "\n",
    "Pour toutes les applications numériques, on choisit : $a=0.1$, $b=d=6\\cdot 10^{-5},$ $c=0.05$. On choisit une population initiale de $H_0=1500$ sardines et $P_0=500$ requins.\n",
    "\n",
    "**EXERCICE 2**\n",
    "\n",
    "> 1. Reprenez chacune des questions des Exercices 1-3 du TP02 pour ce modèle. Visualisez la solution pour les valeurs de p  de la liste $ [0, 0.02, 0.04, 0.06] $. On pourra tracer le champs de vecteurs pour $ p = 0 $ uniquement.\n",
    "> 2. Expliquez quel est l'impact de la pêche sur l'évolution des populations. Ce modèle permet-il d'expliquer les observations de d'Anconna sur la peche de la sardine dans la mer Adriatique (voir le cours)? Justifiez votre réponse.\n",
    "> 3. Quelle espèce est favorisée par la pêche? "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 584664 characters]\n",
      "text/plain": [
       "<Figure size 1200x600 with 3 Axes>"
      ]
     },
     "metadata": {
      "image/png": {
       "height": 547,
       "width": 1014
      }
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Définition des paramètres\n",
    "a, b, c, d = 0.1, 6e-5, 0.05, 6e-5\n",
    "T, N = 200, 100\n",
    "H0, P0 = 1500, 500\n",
    "\n",
    "def F(X, t, a, b, c, d, p):\n",
    "    \"\"\"Fonction second membre pour le système\"\"\"\n",
    "    x, y = X\n",
    "    return np.array([x * (a - p - b*y), y * (-c - p + d*x)])\n",
    "\n",
    "t = np.linspace(0, T, N+1)\n",
    "sardines, requins = np.meshgrid(\n",
    "    np.linspace(0.1, 3000, 20),\n",
    "    np.linspace(0.1, 4500, 30)\n",
    ")\n",
    "fsardines = F((sardines, requins), t, a, b, c, d, 0)[0]\n",
    "frequins = F((sardines, requins), t, a, b, c, d, 0)[1]\n",
    "n_sndmb = np.sqrt(fsardines**2 + frequins**2)  \n",
    "\n",
    "# On crée une figure à trois graphiques\n",
    "fig = plt.figure(figsize=(12, 6))\n",
    "ax = fig.add_subplot(1, 2, 2) # subplot pour le champ de vecteurs et le graphe sardines vs requins\n",
    "axr = fig.add_subplot(2, 2, 1) # subplot pour le graphe du nombre de requins en fonction du temps\n",
    "axs = fig.add_subplot(2, 2, 3) # subplot pour le graphe du nombre de sardines en fonction du temps  \n",
    "ax.quiver(sardines, requins, fsardines/n_sndmb, frequins/n_sndmb)\n",
    "\n",
    "list_p = [0, 0.02, 0.04, 0.06]\n",
    "for k, pk in enumerate(list_p):\n",
    "    couleur = (0, k/len(list_p), 1-k/len(list_p))\n",
    "    X = odeint(F, np.array([H0, P0]), t, args=(a, b, c, d, pk))\n",
    "    \n",
    "    # Tracer la courbe parametrée (H(t),P(t)) \n",
    "    ax.plot(X[:, 0], X[:, 1], linewidth=2, color=couleur, label=f\"$p={pk}$\")\n",
    "    \n",
    "    # Tracer H en fonction du temps \n",
    "    axs.plot(t, X[:, 0], label=f\"Sardines pour p={pk}\", color=couleur)\n",
    "    \n",
    "    # Tracer P en fonction du temps\n",
    "    axr.plot(t, X[:, 1], label=f\"Requins pour p={pk}\", color=couleur)\n",
    "    \n",
    "ax.axis('equal')\n",
    "ax.set_title(\"Champ de vecteur du problème de Lotka-Volterra\")\n",
    "ax.set_xlabel(\"Sardines\")\n",
    "ax.set_ylabel(\"Requins\")\n",
    "ax.legend()\n",
    "\n",
    "axr.set_xlabel(\"Temps t\")\n",
    "axr.set_ylabel(\"Requins\")\n",
    "\n",
    "axs.set_xlabel(\"Temps t\")\n",
    "axs.set_ylabel(\"Sardines\")\n",
    "\n",
    "axs.set_title('Evolution des sardines')\n",
    "axr.set_title('Evolution des requins')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": []
   },
   "source": [
    "## 3. Schéma de Crank-Nicolson\n",
    "\n",
    "Le *schéma de Crank-Nicolson* consiste à construire la suite $ y^n $ des valeurs approchées de $ y( t^n ) $, $ n = 0, 1, \\ldots, N $, via les relations suivantes:\n",
    "\n",
    "\\begin{equation*}\n",
    "\\left\\lbrace\n",
    "\\begin{aligned}\n",
    "& y^{ n+1 } = y^n + \\frac{\\Delta t}{ 2 }\\left( f( t_{ n }, y^{ n } ) + f( t_{ n+1 }, y^{ n+1 } )\n",
    "\\right), \\quad n = 0, 1, 2, \\ldots, N-1, \\\\\n",
    "& y^0 = y_0.\n",
    "\\end{aligned}\n",
    "\\right.\n",
    "\\end{equation*}\n",
    "\n",
    "\n",
    "**EXERCICE 3**\n",
    "> 1. Le schéma de Crank-Nicolson est-il explicite ou implicite? Justifiez votre réponse. \n",
    "> 2. Pour la fonction $ f $ définie par $ f( t, y ) = r y $ écrivez la relation qui donne $y^{n+1}$ en fonction de $y^n$ pour le schéma de Crank-Nicolson.\n",
    "> 3. Programmez la méthode *crank_nicolson* afin de calculer la solution approchée du problème de Cauchy du modèle de Malthus. On écrira une fonction `crank_nicolson` qui prend en argument un nombre `y0` (la valeur initiale de la concentration), un nombre `T` (la longueur de l'intervalle de temps), un nombre $ N $ (le nombre de points) et un nombre `r` (le taux de croissance dans le modèle de Malthus), et qui retourne deux tableaux : `t`, constitué des instants $ 0, \\Delta t, 2 \\Delta t, \\ldots, N \\Delta t $, avec $ \\Delta t = T / N $, et `y`, constitué des valeurs de la solution approchée en ces instants. On étudiera le modèle de Malthus avec $ r = -4 $ et on prendra $ y_0 = 1 $, $ T = 16 $.\n",
    "> 4. Pour la méthode d'Euler explicite et pour le schéma de Crank-Nicolson, tracez sur un même graphique, la solution exacte et les solutions approchées du modèle de Malthus pour les valeurs $ N = 31, 32, 33, 34$. \n",
    "Commentez le graphique obtenu. Pour trouver une solution approchée du modèle de Malthus étudié quelle méthode numérique choisiriez-vous ? Expliquez votre choix. \n",
    "> 5. Pour l'équation de Malthus, calculez et affichez les erreurs globales du schéma de Crank-Nicolson pour différentes valeurs de $\\Delta t$. Estimez l'ordre de convergence du schéma de Crank-Nicolson pour l'équation de Malthus. L'ordre de convergence du schéma de Crank-Nicolson est-il plus elevé que l'ordre de convergence du schéma d'Euler implicite?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "def crank_nicolson(y0, T, N, r):\n",
    "    \"\"\"\n",
    "    schéma de Crank-Nicolson pour le modèle de Malthus \n",
    "    \n",
    "    Parameters\n",
    "    ----------\n",
    "    \n",
    "    y0: float\n",
    "        donnée initiale\n",
    "    T: float\n",
    "        instant final\n",
    "    N: int\n",
    "        nombre de pas de temps\n",
    "    r: float\n",
    "       taux de croissance dans le modèle de Malthus\n",
    "\n",
    "    Returns\n",
    "    -------\n",
    "    \n",
    "    t: ndarray\n",
    "        les instants où la solution approchée est calculée\n",
    "    y: ndarray\n",
    "        les valeurs de la solution approchée par le theta-schema\n",
    "    \"\"\"\n",
    "    \n",
    "    dt = T / N\n",
    "    t = np.zeros(N+1)\n",
    "    y = np.zeros(N+1)\n",
    "    tk, yk = 0, y0\n",
    "    y[0] = yk\n",
    "    \n",
    "    for n in range(N):\n",
    "        tk += dt\n",
    "        yk *= (2 + dt * r) / (2 - dt * r)\n",
    "        y[n+1] = yk\n",
    "        t[n+1] = tk\n",
    "        \n",
    "    return t, y\n",
    "\n",
    "def euler_explicit(y0, T, N, r):\n",
    "    \"\"\"\n",
    "    schéma de d'Euler pour le modèle de Malthus \n",
    "    \n",
    "    Parameters\n",
    "    ----------\n",
    "    \n",
    "    y0: float\n",
    "        donnée initiale\n",
    "    T: float\n",
    "        instant final\n",
    "    N: int\n",
    "        nombre de pas de temps\n",
    "    r: float\n",
    "       taux de croissance dans le modèle de Malthus\n",
    "\n",
    "    Returns\n",
    "    -------\n",
    "    \n",
    "    t: ndarray\n",
    "        les instants où la solution approchée est calculée\n",
    "    y: ndarray\n",
    "        les valeurs de la solution approchée par le theta-schema\n",
    "    \"\"\"\n",
    "    dt = T / N\n",
    "    t = np.zeros(N+1)\n",
    "    y = np.zeros(N+1)\n",
    "    tk, yk = 0, y0\n",
    "    y[0] = yk\n",
    "    \n",
    "    for n in range(N):\n",
    "        tk += dt\n",
    "        yk += dt * r * yk\n",
    "        y[n+1] = yk\n",
    "        t[n+1] = tk\n",
    "        \n",
    "    return t, y\n",
    "\n",
    "def solution_exacte(t):\n",
    "    \"\"\"\n",
    "    Fonction calculant la solution exacte du modèle de Malthus à l'instant t \n",
    "    \"\"\"\n",
    "    return y0 * np.exp(r * t)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0, 0.5, 'y')"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "[encoded data, 202120 characters]\n",
      "text/plain": [
       "<Figure size 1600x600 with 2 Axes>"
      ]
     },
     "metadata": {
      "image/png": {
       "height": 547,
       "width": 1309
      }
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Définition des paramètres\n",
    "r = -4\n",
    "liste_N = [31, 32, 33, 34]\n",
    "T = 16\n",
    "y0 = 1\n",
    "\n",
    "t_exact = np.linspace(0, T, 1000)\n",
    "\n",
    "fig = plt.figure(figsize=(16, 6))\n",
    "\n",
    "# Schéma d'Euler explicite\n",
    "ax = fig.add_subplot(1, 2, 1)\n",
    "for n in liste_N:\n",
    "    t, y = euler_explicit(y0, T, n, r) # On calcule la fonction Euler pour chaque n\n",
    "    ax.scatter(t, y, label=f\"Solution approchée pour N={n}\")\n",
    "    \n",
    "ax.plot(t_exact, solution_exacte(t_exact), label='Solution exacte')\n",
    "ax.legend()\n",
    "ax.axis('equal')\n",
    "ax.set_title(\"Schéma d'Euler explicite\")\n",
    "ax.set_xlabel(\"Temps t\")\n",
    "ax.set_ylabel(\"y\")\n",
    "    \n",
    "\n",
    "# Schéma de Crank-Nicolson\n",
    "ax = fig.add_subplot(1, 2, 2)\n",
    "for n in liste_N:\n",
    "    t, y = crank_nicolson(y0, T, n, r) # On calcule la fonction Crank-Nicolson pour chaque n\n",
    "    ax.scatter(t, y, label=f\"Solution approchée pour N={n}\")\n",
    "ax.plot(t_exact, solution_exacte(t_exact), label='Solution exacte')\n",
    "ax.legend()\n",
    "ax.set_title(\"Schéma de Crank-Nicolson\")\n",
    "ax.set_xlabel(\"Temps t\")\n",
    "ax.set_ylabel(\"y\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Delta_t =  1.600e-01, e =  5.255e-01\n",
      "Delta_t =  1.600e-01, e =  2.715e-01\n",
      "Delta_t =  1.600e-01, e =  5.853e-02\n",
      "Delta_t =  1.600e-01, e =  1.266e-02\n",
      "Delta_t =  1.600e-01, e =  3.173e-03\n",
      "Delta_t =  1.600e-01, e =  5.030e-04\n",
      "Delta_t =  1.600e-01, e =  1.256e-04\n"
     ]
    },
    {
     "data": {
      "image/png": "[encoded data, 100068 characters]\n",
      "text/plain": [
       "<Figure size 600x600 with 1 Axes>"
      ]
     },
     "metadata": {
      "image/png": {
       "height": 547,
       "width": 547
      }
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "liste_N = np.array([10, 20, 50, 100, 200, 500, 1000])\n",
    "liste_delta = T / liste_N\n",
    "liste_erreur = np.zeros(liste_N.size)\n",
    "\n",
    "for k, n in enumerate(liste_N):\n",
    "    t, sol_appr = crank_nicolson(y0, T, n, r)\n",
    "    sol_ex = solution_exacte(t)\n",
    "    erreur = np.max(np.abs(sol_appr - sol_ex))\n",
    "    #erreur = np.linalg.norm(sol_appr - sol_ex, np.inf)\n",
    "    print(f\"Delta_t = {T / N:10.3e}, e = {erreur:10.3e}\")\n",
    "    liste_erreur[k] = erreur\n",
    "\n",
    "fig = plt.figure(figsize=(6, 6))\n",
    "ax = fig.add_subplot(1, 1, 1)\n",
    "\n",
    "ax.scatter(liste_delta, liste_erreur, color=\"black\")\n",
    "for p in [0.5, 1, 2]:\n",
    "    C = liste_erreur[-1] / (liste_delta[-1] ** p)\n",
    "    plt.plot(liste_delta, C * liste_delta ** p, label=f\"$p={p}$\")\n",
    "ax.set_title(\"Erreur du schéma de Crank-Nicolson\")\n",
    "ax.set_xlabel(r\"$\\Delta t$\")\n",
    "ax.set_ylabel(r\"$e(\\Delta t)$\")\n",
    "ax.set_xscale(\"log\")\n",
    "ax.set_yscale(\"log\")\n",
    "ax.legend()\n",
    "plt.show()"
   ]
  },
  {
   "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"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
