{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "d54b4116-d2ec-4cbf-b914-68a322028a25",
   "metadata": {},
   "source": [
    "# Résolution numérique d’équations différentielles"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "308b01b9",
   "metadata": {},
   "source": [
    "> **En fin de séance, rendre le fichier source Jupiter Notebook (`.ipynb`) de votre TP obtenu <span style=\"color:red\">après les 3 heures</span> de travail.**\n",
    "> \n",
    ">**Rendre précisément <span style=\"color:red\">une semaines après</span> la séance de TP le fichier source final comportant les réponses aux questions (dans ce TP, vous avez trois réponses à rédiger pour chaque exercice) et les calculs demandés.**\n",
    "> \n",
    "> L'objectif de la seance est d'étudier des **méthodes numériques** permettant obtenir des solutions d’équations différentielles.\n",
    ">\n",
    "> Vous allez comparer les solutions théoriques obtenues à l’aide des notions vu en CM avec les solutions qu’on peut obtenir à l'aide de ces méthodes.\n",
    ">\n",
    "> Par ailleurs vous allez apprendre la redaction des formules mathématiques en $\\mathrm\\LaTeX$ pour écrire vos réponses aux questions ci-dessous (un petit manuel est disponible sur moodle pour vous aider).\n",
    "\n",
    "**Ces travaux, contenant vos réponses aux questions, seront utilisés pour évaluer les compétences suivantes :**\n",
    "- Identifier un modèle mathématique lors de la mise en équation d'une expérience scientifique.\n",
    "- Effectuer un calcul élémentaire à l'aide d'un logiciel de calcul symbolique ou numérique.\n",
    "- Interpréter les résultats obtenus afin de répondre à une problématique.\n",
    "- Utiliser les outils de communication courants.\n",
    "- Structurer/rédiger/présenter selon une méthodologie adaptée.\n",
    "    \n",
    "**Cette évaluation comptera pour 20 % de votre note finale dans ce cours.**"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1945a7aa",
   "metadata": {},
   "source": [
    "## Exercice 1. Modèle de croissance (équations différentielles appliquées à la biologie)\n",
    "\n",
    "Nous allons d'abord étudier une équation du premier ordre et **la méthode d’Euler** pour construire la solution numérique. \n",
    "\n",
    "Il existe de nombreux modèles de croissance. Considérons par exemple le modèle suivant :\n",
    "$$ N' = (2 - \\alpha\\cos(t))N - F(N) \\tag{1} $$\n",
    "où\n",
    "- $N(t)$ est l’effectif de la population en fonction du temps $t$ ;\n",
    "- $2 - \\alpha\\cos(t)$ est le taux de naissance avec variations saisonnières, le coefficient $\\alpha$ (alpha) est constant ;\n",
    "- $F(N) = \\frac{\\beta}{2} N^2$ est le terme de surpopulation qui depend de la solution elle-même, le coefficient $\\beta$ (beta) est constant.\n",
    "\n",
    "Comme le terme $F(N)$ est une fonction de $N$, cette équation est dite non linéaire. En général, on ne peut pas toujours construire explicitement une solution à une équation non linéaire. D’où l'utilité des méthodes numériques.\n",
    "\n",
    "Dans le cas $F(N) = 0$, on obtient une équation homogène à coefficients variables, contrairement aux exemples étudiés en CM. Cela rend plus difficile la résolution du problème.\n",
    "\n",
    "### Travail à realiser\n",
    "\n",
    "1. Montrer que, dans le cas homogène $F(N) = 0$, la solution de l’équation (1) est donnée par la relation\n",
    "$$ y(t) = K \\exp(2t - \\alpha \\sin t), \\quad K = \\text{const}. $$\n",
    "Vous expliciterez vos calculs dans la cellule ci-dessous.\n",
    " \n",
    "    (Rédigez votre reponse en $\\mathrm\\LaTeX$. **Cette réponse est obligatoire pour le rendu dans une semaine.**)\n",
    "\n",
    "    *Indication. Réécrire l’équation sous la forme\n",
    "    $$ \\frac{y'}{y}=a(t) $$\n",
    "    et intégrer les deux cotés. On admettra que la fonction $y$ ne s’annule jamais.*"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d0c3eb4b",
   "metadata": {},
   "source": [
    "> 1. **Votre réponse ici. (Double-cliquer ici pour modifier le contenu.)**  \n",
    "Pour verifier que la fonction $y(x)$ est une solution, on fait le calcul..."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8cf3cd4c",
   "metadata": {},
   "source": [
    "2. Étudiez l’exemple d’utilisation de la méthode `plt.plot` pour tracer les fonctions. Dans l’exemple, le nombre de points dans le tableau `x` est imposé à 1000. Quel est le nombre minimal pour avoir un graphique toujours bien presenté ? Vous changerez la valeur dans l'exemple."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d9e90a99",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# En python, on commence toujours par importer les bibliothèques nécessaires\n",
    "# la bibliothèque matplotlib.pyplot est utilisée pour tracer des fonctions\n",
    "# la bibliothèque numpy permet de travailler efficacement avec des tableaux \n",
    "#                           et est souvent utilisée pour les calculs scientifiques.\n",
    "\n",
    "\n",
    "import matplotlib.pyplot as plt \n",
    "import numpy as np\n",
    "import math as math\n",
    "\n",
    "# les instructions ci-dessous permettent de definir des focntions f1 et f2 d'un argument x\n",
    "# l'utilisation de np.sin et np.exp (méthodes de la bibliothèque numpy) permet\n",
    "#                                  d'utiliser ces fonctions même dans le cas où x est un tableau\n",
    "\n",
    "def f1(x):\n",
    "    return x * np.sin(3 * x**2)\n",
    "\n",
    "def f2(x):\n",
    "    return np.exp(-x**2 / 2)\n",
    "    \n",
    "# la commande np.linspace permet de créer un tableau de 1000 points equirepartis entre -pi et pi\n",
    "x = np.linspace(-math.pi, math.pi, 1000)\n",
    "# la ligne ci-dessous permet de créer un tableau y avec des valeurs de la fonction f1\n",
    "y = f1(x)\n",
    "\n",
    "# méthode plt.plot() pour tracer les points (x, y)\n",
    "plt.plot(x, y)\n",
    "# plt.show() pour afficher le graphique\n",
    "plt.show()\n",
    "\n",
    "# on peut utiliser f1 et f2 directiment quand on appele la méthode plt.plot\n",
    "plt.plot(x, f1(x), 'r')\n",
    "plt.plot(x, f2(x), 'b')\n",
    "# si maintenant on appelle plt.show() juste après, on peut afficher les deux graphiques ensemble\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5a21e2e9",
   "metadata": {},
   "source": [
    "Définissez une fonction périodique $f_3$ et représentez-la graphiquement sur deux périodes."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6fd7f1e3",
   "metadata": {},
   "outputs": [],
   "source": [
    "def f3(x):\n",
    "    return "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "96350b4f",
   "metadata": {},
   "source": [
    "3. Revenons à notre équation : modifiez l’exemple ci-dessus pour représenter graphiquement une solution de l’équation homogène de la question 1."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "47af2450",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Graphique d'une solution de l'equation homogène \n",
    "\n",
    "def yh(t):\n",
    "    return\n",
    "\n",
    "# t = ..\n",
    "# y = ...\n",
    "# plt.plot(...)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b9e43c80",
   "metadata": {},
   "source": [
    "4. On passe maintenant à la résolution numérique de l’équation (1). Pour cela on va utiliser **la méthode d’Euler**.\n",
    "La méthode d’Euler permet d'obtenir une solution approchée de l’équation différentielle\n",
    "$$\n",
    "\\left\\{ \\begin{aligned}\n",
    "    y' &= f(x, y), \\\\\n",
    "    y(x_0) &= y_0.\n",
    "\\end{aligned} \\right.\n",
    "$$\n",
    "Il est à noter que la donnée $y_0$ est connue. Ici, la lettre $f$ désigne une fonction de deux variables à valeurs réelles, c'est-à-dire\n",
    "$$\n",
    "f \\colon \\left\\vert \\begin{aligned}\n",
    "    \\mathbb R \\times \\mathbb R &\\longrightarrow \\mathbb R, \\\\\n",
    "    (x, y) &\\longmapsto f(x, y).\n",
    "\\end{aligned} \\right.\n",
    "$$\n",
    "Pour l'équation (1), on a\n",
    "$$ f(t, N(t)) = (2 - \\alpha\\cos(t))N(t) - \\frac{\\beta}{2}N(t)^2. $$\n",
    "On fixe un pas d'itération $h$ et on définit la suite $(t_0, \\dots, t_M)$ par la relation\n",
    "$$\n",
    "t_j = t_0 + jh, \\quad j = 0, \\dots, M.\n",
    "$$\n",
    "Les valeurs approchées $N_j$ de $N(t_j)$ sont définies itérativement en posant\n",
    "$$\n",
    "N_{j + 1} = N_j + hf(t_j, N_j). \\tag{2}\n",
    "$$\n",
    "Complétez l’algortithme Python ci-dessous pour le calcul des $M$ itérations de la suite $(N_0, \\dots, N_M)$ définie à l’aide de la méthode d’Euler avec un pas $h$ et une donnée initiale $N_0$ (correspondant à $N(t_0)$ avec $t_0 = 0$) qu’on pourrait modifier à notre convenance. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "91404668",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Méthode d’Euler pour la résolution de l’équation (1)\n",
    "\n",
    "def Euler(h, N0, M, alp, beta):\n",
    "    N = np.zeros(M)\n",
    "    t = np.zeros(M)\n",
    "    N[0] = ...\n",
    "    t[0] = 0\n",
    "    for j in ...\n",
    "        N[j + 1] = N[j] + h...\n",
    "        t[j + 1] = t[j] + h...\n",
    "    return N, t"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "91147142",
   "metadata": {},
   "source": [
    "5. Comparez graphiquement les valeurs correspondantes à la solution exacte obtenue à la question 1 (choisissez vos paramètres avec soin !) et les valeurs calculées numériquement. Changez les valeurs du pas $h$ et étudier son influence sur les résultats."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d86f42d4",
   "metadata": {},
   "source": [
    "6. Comment le choix du pas $h$ affecte-t-il les résultats ? Commentez."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2d9232b0",
   "metadata": {},
   "source": [
    "> 6. **Votre réponse ici. (Double-cliquer ici pour modifier le contenu)**  <br>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "08a2b6f2",
   "metadata": {},
   "source": [
    "7. Que pouvez vous dire de comportement de la solution dans le cas homogène ? Que se passe-t-il avec la population ?"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5eb7fe28",
   "metadata": {},
   "source": [
    "> 7. **Votre réponse ici. (Double-cliquer ici pour modifier le contenu)**"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d295c163",
   "metadata": {},
   "source": [
    "8. Construire maintenant des solution pour les differentes cas où $\\beta \\neq 0$ et representer les graphiquement. Commentez le résultat. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "278e6123",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "0c3cf0bd",
   "metadata": {},
   "source": [
    "### Exercice 2. Régimes d'évolution d'oscillateurs libres amorties\n",
    "\n",
    "On sait dèja que le mouvement d'un pendule se decrit par une équation de second ordre. On considère l’équation suivante qui permet de prendre en compte l'amortissement :\n",
    "$$\n",
    "y''(t) + 2\\gamma y'(t) + \\omega_0^2y(t) = 0 \\tag{3}\n",
    "$$\n",
    "où $\\gamma > 0$ correspond à l’amortissement et $\\omega_0 > 0$ correspond à la fréquence propre.\n",
    "\n",
    "1. Donnez la forme des solutions de l’équation différentielle ci-dessus. On distinguera trois cas de figure selon la valeur du nombre $\\Delta = 4(\\gamma^2 - \\omega_0^2)$. (Vous pouvez revenir sur l'exercice 2 de TD 1, seulement ici on va considerer les solutions réelles directement.)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f34ee9db",
   "metadata": {},
   "source": [
    "> 1. **Votre réponse ici. (Rédigez votre reponse en $\\mathbf\\LaTeX$.)**"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a3c384c8",
   "metadata": {},
   "source": [
    "2. La méthode d'Euler peut également être utilisée pour résoudre des équations du second ordre si l'on peut les réduire à un système de deux équations du premier ordre auxquelles on applique la méthode comme dans l'exercice précédent, notamment pour une équation avec les données $y_0$ et $z_0$ comme\n",
    "$$\n",
    "\\begin{aligned}\n",
    "    y''(t) + f(t, y(t), y'(t)) &= 0, \\\\\n",
    "    y(0) &= y_0, \\\\\n",
    "    y'(0) &= z_0.\n",
    "\\end{aligned}\n",
    "$$\n",
    "On introduit la variable $z(t) = y'(t)$ et on obtient un système\n",
    "$$\n",
    "\\left\\{ \\begin{aligned}\n",
    "    y'(t) &= z(t), \\\\\n",
    "    z'(t) &= -f(t, y(t), z(t)), \\\\\n",
    "    y(0) &= y_0, \\\\\n",
    "    z(0) &= z_0.\n",
    "\\end{aligned} \\right.\n",
    "$$\n",
    "On fixe encore un pas d'itération $h$ et on définit la suite $(t_0, \\dots, t_M)$ par la relation\n",
    "$$\n",
    "t_j = t_0+j h, \\quad j= 0\\ldots M.\n",
    "$$\n",
    "où la quantité $t_0$ est une donnée. Dans notre cas, on prend $t_0 = 0$. La méthode d'Euler s’écrit \n",
    "$$\n",
    "\\left\\{ \\begin{aligned}\n",
    "    y_{j + 1} &= y_j + hz_j, \\\\\n",
    "    z_{j + 1} &= z_j - hf(t_j, y_j, z_j),\n",
    "\\end{aligned} \\right.\n",
    "$$\n",
    "où on itialise les suites $(y_0, \\dots, y_M)$ et $(z_0, \\dots z_M)$ respectivement par les données $y_0$ et $z_0$.\n",
    "\n",
    "    Écrivez explicitement la méthode d’Euler appliquée à l'équation (3) : mettez d’abord l'équation (3) sous la forme d'un système et explicitez ensuite la méthode pour les itérés des valeurs approchées de $y(x_j)$ et $y'(x_j)$."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "be1206ff",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "source": [
    "> 2. **Votre réponse ici. (Rédigez votre reponse en $\\mathbf\\LaTeX$.)**"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9f2b9755",
   "metadata": {},
   "source": [
    "3. Implémentez l’algorithme Python pour la methode d'Euler pour la résolution de l’équation (3) en adoptant votre méthode de la question 4 de l'exercice précédent."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "625a8414",
   "metadata": {},
   "outputs": [],
   "source": [
    "def EulerP(h, y0, yp0, M, g, w0):\n",
    "    t = np.zeros(M)\n",
    "    y = np.zeros(M)\n",
    "    z = np.zeros(M)\n",
    "    t[0] = 0\n",
    "    y[0] = y0\n",
    "    z[0] = yp0\n",
    "    for j in range(M - 1):\n",
    "        z[j + 1] = ...\n",
    "        y[j + 1] = ...\n",
    "        t[j + 1] = ...\n",
    "    return y, x"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2fb97351",
   "metadata": {},
   "source": [
    "On distingue quatre régimes possibles pour ce système physique :\n",
    "\n",
    "- régime apériodique : $\\Delta > 0$. On dit que ce sont des solutions amorties ou surcritiques, il n’y a pas d’oscillation autour de l’axe des abscisses ;  \n",
    "- régime apériodique critique : $\\Delta = 0$. On dit que ce sont des solutions amorties critiques, il n’y a pas d’oscillation autour de l’axe des abscisses ;  \n",
    "- régime pseudo-périodique ou pseudo-critique : $\\Delta < 0$ et $\\gamma \\neq 0$. Il y a donc des oscillations autour de l’axe des abscisses de moins en moins grande ;  \n",
    "- régime harmonique : $\\Delta < 0$ et $\\gamma = 0$. Il s’agit d'un régime théorique (exemple considéré dans le CM) car normalement, dans tout système physique, il y a des pertes d'énergie (frottement)."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f5f516c7",
   "metadata": {},
   "source": [
    "4. Choisissez bien les valeurs des paramètres et représenter le déplacement $y(x)$ obtenu numériquement à l’aide de votre algorithme d’Euler pour tous les régimes possibles.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "af5bd197",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "6a5eab6c",
   "metadata": {},
   "source": [
    "5. On se concetre maintenant sur le cas de régime harmonique $\\gamma = 0$. Représentez sur le même graphique la solution exacte et la solution obtenue à l’aide de votre méthode numérique. Faites varier la valeur de $h$ d'un petit pas à un grand pas. Commentez les résultats observés. Vous allez observer des problèmes liés aux propriétés de la méthode d’Euler (et, en particulier, à sa stabilité)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4705777a",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "dd21fa48",
   "metadata": {},
   "source": [
    "6. Peut-on améliorer la méthode ? Il est possible d’utiliser une amélioration de la méthode d’Euler, appellée la méthode de Crank-Nicolson, qui se presente de la façon suivante :\n",
    "$$\n",
    "\\left\\{ \\begin{aligned}\n",
    "    y_{j + 1} &= y_j + h\\frac{z_j + z_{j + 1}}{2}, \\\\\n",
    "    z_{j + 1} &= z_j - h\\frac{f(x_j, y_j, z_j) + f(x_{j + 1}, y_{j+1}, z_{j + 1})}{2}.\n",
    "\\end{aligned} \\right.\n",
    "$$\n",
    "Selon la fonction $f$, nous pouvons ou non déterminer les valeurs $y_{j + 1}$ et $z_{j + 1}$ explicitement. Dans notre cas, c’est possible. Écrivez explicitement la méthode de Crank-Nicolson appliquée à l'équation (3) : mettez d'abord l'équation (3) sous la forme d'un système comme précédemment, puis exprimez les itérés des valeurs approchées $y_j$ et $z_j$ de $y(x_j)$ et $y'(x_j)$."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9653cdee",
   "metadata": {},
   "source": [
    "> 6. **Votre réponse ici. (Rédigez votre reponse en $\\mathbf\\LaTeX$.)**"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b99876fa",
   "metadata": {},
   "source": [
    "7. Implémentez l’algorithme Python pour la résolution de l’équation (3) à l’aide de la méthode de Crank-Nicolson. Refaites les tests pour les mêmes valeurs du pas $h$ pour lesquelles vous avez observé les problèmes dans la question 5."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0ac86bae",
   "metadata": {},
   "outputs": [],
   "source": [
    "def CrankNicolson(h, y0, yp0, M, g, w0):\n",
    "    x = np.zeros(M)\n",
    "    y = np.zeros(M)\n",
    "    z = np.zeros(M)\n",
    "    x[0] = 0\n",
    "    y[0] = y0\n",
    "    z[0] = yp0\n",
    "    for j in range(M - 1):\n",
    "        z[j + 1] = ...\n",
    "        y[j + 1] = ...\n",
    "        x[j + 1] = ...\n",
    "    return y, x"
   ]
  }
 ],
 "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.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}