{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "ef674ab4-d343-474a-9b6e-a3888c353796",
   "metadata": {},
   "source": [
    "## MT09 - TP7 - Automne 2024\n",
    "### Schémas numériques pour les équations différentielles\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3fe9daf9-8d95-4546-9983-987fbaaa6b56",
   "metadata": {
    "tags": []
   },
   "source": [
    "### 1. Croissance cellulaire\n",
    "\n",
    "Dans les exercices de travaux dirigés, on a vu le modèle différentiel continu de croissance cellulaire\n",
    "$$\n",
    "\\frac{d\\rho}{dt} = \\sigma \\frac{\\rho}{\\rho_M}(\\rho_M-\\rho),\\quad t\\in [0,T]\n",
    "$$\n",
    "\n",
    "où $\\rho$ représente la densité de cellules par unité de surface. À cette équation différentielle, on associe la condition initiale $\\rho(0) = \\rho_0$, $0<\\rho_0 \\ll \\rho_M$. On rappelle que la solution exacte s'écrit\n",
    "\n",
    "$$\n",
    "\\rho(t) = \\rho_M \\, \\frac{e^{\\sigma t}}{\\frac{\\rho_M}{\\rho_0}-1+e^{\\sigma t}}.\n",
    "$$\n",
    "\n",
    "\n",
    "Le module ```scipy.integrate``` donne accès à la méthode ```odeint()``` qui résout numériquement un problème différentiel :\n",
    "\n",
    "```\n",
    "y = odeint(model, y0, t)\n",
    "```\n",
    "\n",
    "(voir l'aide de ```odeint()```). Ecrire une fonction\n",
    "```croissanceCell(y,t)``` qui implémente la fonction de l'équation différentielle de croissance cellulaire. Résoudre ensuite numériquement le problème différentiel et comparer à la solution exacte $\\rho(t)$.\n",
    "\n",
    "Pour l'application numérique, on utilisera les paramètres\n",
    "$T=20$, $\\rho_M=1000$, $\\rho_0=5$ et $\\sigma=0.7$ et ```t = np.linspace(0,T, 100)```."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "526623a5-7281-4d52-a273-b2a21957284c",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "def rho(t):\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b3abacb2-af35-4f4a-b671-9b6f3678bd8e",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "from scipy.integrate import odeint\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "# function that returns ydot\n",
    "def croissanceCell(y,t):\n",
    "    sigma = 0.7\n",
    "    rhoM = 1000\n",
    "    ydot = #...\n",
    "    return ydot\n",
    "\n",
    "# initial condition\n",
    "y0 = np.array([5.0])\n",
    "# final time\n",
    "T = 20.0\n",
    "# time points\n",
    "tref = np.linspace(0,T, 100)\n",
    "# solve ODE\n",
    "yref = odeint(croissanceCell, y0, tref)\n",
    "\n",
    "# plot results\n",
    "#..."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "039b31c2-1318-4f52-8a48-380ef43cd036",
   "metadata": {},
   "source": [
    "#### Schéma d'Euler explicite\n",
    "Ecrire une fonction python ```Euler(F, y0, t0, T, N)``` qui met en oeuvre le schéma d'Euler explicite générique pour une fonction $F$ (à valeurs scalaires ou vectorielles), une donnée initiale $y_0$, un temps initial $t_0$, un temps final $T$, sur une subdivision uniforme constituée de $N$ sous-intervalles, avec un pas\n",
    "$h=\\frac{T-t_0}{N}$. Les tableaux de sorties seront le tableau des temps discrets $t$ et le tableau (ou matrice) $z$ des états du système $y_n$ aux temps discrets $t_n$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9eaabc75-4fbd-473e-91f3-9102b2733031",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "def Euler(F, y0, t0, T, N):\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e5a14aa1-5d27-4dc3-be95-799d28da8de7",
   "metadata": {},
   "source": [
    "Appliquez la méthode d'Euler au problème de croissance cellulaire. On prendra toujours $y_0=5$, $T=20$, et $N=50$. Tracer la solution numérique $t\\mapsto z(t)$ et comparer à la solution de référence obtenue avec ```odeint()```. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "320268fa-31c1-4102-890c-b87006abf763",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "y0 = np.array([5.0])\n",
    "T = 20.0\n",
    "N = 50\n",
    "t, z = Euler(croissanceCell, y0, 0.0, T, N)\n",
    "#..."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b2edd6b6-9894-4c1e-94bf-6fbfc712e1c0",
   "metadata": {},
   "source": [
    "Testez plusieurs valeurs de $N$, par exemple $N=10,\\ 50,\\ 100, 200$ pour observer l'évolution de l'erreur en fonction de $N$ et vérifier la convergence de la méthode."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ba825e2e-e89f-456b-8e62-395dc4b9eeb9",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "0ba9aa89-764c-44f0-b7a7-2d6084fc76bc",
   "metadata": {},
   "source": [
    "### 2. Oscillateur non linéaire de Duffing\n",
    "\n",
    "Le système non linéaire de Duffing modélise un système masse-ressort isolé avec une loi d'élasticité non linéaire cubique. Les équations de la cinématique et de la dynamique autour de l'équilibre donnent\n",
    "\n",
    "\\begin{align*}\n",
    "& \\dot x(t) = v(t), \\\\\n",
    "& \\dot v(t) = - \\alpha\\, x(t) -\\beta\\, x^3(t).\n",
    "\\end{align*}\n",
    "\n",
    "Le terme $\\alpha x$ est le terme d'élasticité linéaire classique avec $\\alpha>0$ le coefficient de raideur du ressort. Le terme non linéaire\n",
    "$\\beta x^3$, $\\beta>0$, agit pour les plus forts déplacements et exprime une plus grande force de rappel.\n",
    "\n",
    "On peut vérifier que l'énergie totale du système (somme d el'énergie élastique et de l'énergie cinétique)\n",
    "\n",
    "$$\n",
    "E=E(x,v) = \\alpha \\frac{x^2}{2} + \\beta \\frac{x^4}{4}\n",
    "+ \\frac{v^2}{2}\n",
    "$$\n",
    "\n",
    "est conservée au cours du temps (système hamiltonien).\n",
    "\n",
    "Ecrire d'abord le système sous la forme vectorielle générique\n",
    "\n",
    "$$\n",
    "\\frac{d\\mathbf{y}}{dt} = \\mathbf{F}(t, \\mathbf{y}(t))\n",
    "$$\n",
    "\n",
    "où l'on précisera $\\mathbf{y}$ et $\\mathbf{F}(t,\\mathbf{y})$. Ecrire une fonction python ```Duffing(y,t)``` qui implémente $F$.\n",
    "Pour les paramètres, on prendra $\\alpha=1$ et $\\beta=12$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5f77a176-60cb-4158-bbc7-85ac797be179",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# function that returns ydot\n",
    "def Duffing(y,t):\n",
    "\n",
    "    return ydot"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "eb5c0628-5ff4-4d94-9517-7d480384ee14",
   "metadata": {
    "tags": []
   },
   "source": [
    "Intégrer numériquement le problème différentiel en utilisant ```odeint()```. On prendra $t_0=0$, $T=4$, \n",
    "$y_0=(x_0,v_0)^T=(2,0)^T$ et ```t = np.linspace(0,T, 200)```.\n",
    "Tracer les solutions numériques $t\\mapsto x(t)$ et $t\\mapsto v(t)$ obtenues. Tracer de même l'évolution de l'énergie totale $t\\mapsto E(t)$ et faites attention à l'échelle de l'axe des ordonnées."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d82c7fb8-e7ca-4947-b419-71240bcf2a32",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from scipy.integrate import odeint\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "# initial condition\n",
    "y0 = np.array([2, 0])\n",
    "# final time\n",
    "T = 4\n",
    "# time points\n",
    "t = np.linspace(0,T, 200)\n",
    "# solve ODE\n",
    "y = odeint(Duffing, y0, t)\n",
    "#..."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7d3b9136-7065-4677-b97b-9143f4600960",
   "metadata": {},
   "source": [
    "Tracer de même la solution $(x(t),v(t))_{t\\in[0,T]}$ dans le plan de phase $(x,v)$. La courbe fermée ainsi obtenue est dite une orbite."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3e3081a6-0eb7-4b1b-8fc5-95cee3abb261",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "7a4e2a84-b2fd-4886-9d98-08ece9f06110",
   "metadata": {},
   "source": [
    "Appliquer le schéma d'Euler à l'oscillateur de Duffing avec $N=500$. Tracer les solutions numériques obtenues ainsi que l'évolution de l'énergie totale. Est-ce satisfaisant ? "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c49bf475-54a7-42ee-85a8-ccd109961c8e",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8d64568e-bf66-4400-9277-2213b3903796",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "1ec97e7b-bd9d-42e3-b7b8-3f3585c2ceea",
   "metadata": {},
   "source": [
    "#### Schéma de Runge-Kutta RK2\n",
    "\n",
    "Le schéma de Runge-Kutta RK2 (ou schéma d'Euler-Cauchy) du cours s'écrit\n",
    "\n",
    "\\begin{align*}\n",
    "& k_1 = F(t_n, z_n), \\\\[1.3ex]\n",
    "& k_2 = F(t_n+h, z_n+h\\,k_1) \\\\\n",
    "& z_{n+1} = z_n + \\frac{h}{2}\\left( k_1+k_2\\right).\n",
    "\\end{align*}\n",
    "\n",
    "Mettre en oeuvre ce schéma dans une méthode ```RK2(F, y0, t0, T, N)```.\n",
    "Appliquez RK2 au système de Duffing où l'on comparerea les solutions avec celles obtenues avec ```odeint()``` ainsi que l'évolution de l'énergie totale. Commentez. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "622dbff9-c7a7-42b0-a2d2-9276bd3b91fa",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# (ACU)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0323a455-b410-4d74-8483-3e1573f43648",
   "metadata": {},
   "source": [
    "#### Schéma à variables décalées de Verlet\n",
    "\n",
    "\\begin{align*}\n",
    "\t& x_{n+1} = x_n + h v_{n+1/2}, \\\\\n",
    "\t& v_{n+3/2} = v_{n+1/2} - h\\alpha\\, x_{n+1} - h \\beta\\, x_{n+1}^3.  \n",
    "\\end{align*}\n",
    "\n",
    "Les données initiales sont toujours $x_0=2$ et $v_0=0$. Pour \"l'initialisation\" de $v_{1/2}$ on utilisera ici simplement le schéma d'Euler explicite sur $v$ (sachant que l'on peut faire mieux) :\n",
    "\n",
    "$$\n",
    "v_{1/2} = v_0 - \\frac{h}{2}\\alpha\\, x_0 - \\frac{h}{2}\\beta\\, x_0^3.\n",
    "$$\n",
    "\n",
    "Mettre en oeuvre le schéma de Verlet pour l'oscillateur de Duffing. On utilisera une subdivision uniforme comme pour le schéma d'Euler avec $h=\\frac{T-t_0}{N}$. Tracer les solutions numériques $t\\mapsto x(t)$ et $t\\mapsto v(t)$ ainsi obtenues. On utilisera \n",
    "$N=100$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4409d7d5-e980-40db-8004-962b7c2bbaf0",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "3d4f9f88-7a43-4812-97d7-f2a6eb876be1",
   "metadata": {},
   "source": [
    "Les variables discrètes de positions $x_n$ et de vitesse $v_{n+1/2}$ ne sont pas définies aux mêmes instants discrets, si bien que l'on ne peut pas définir strictement une énergie totale. On peut toutefois définir une estimation en calculant une valeur moyenne $\\hat v_n$ de vitesse aux temps $t_n$ par\n",
    "\n",
    "$$\n",
    "\\hat v_n = \\frac{v_{n-1/2}+v_{n+1/2}}{2}\n",
    "$$\n",
    "\n",
    "et définir une énergie totale approchée $\\hat E_n$ par\n",
    "\n",
    "$$\n",
    "\\hat E_n = \\alpha\\frac{x_n^2}{2}+\\beta \\frac{x_n^4}{4} + \\frac{\\hat v_n^2}{2}.\n",
    "$$\n",
    "\n",
    "Tracer dans ce cas l'évolution de $\\hat E_n$ ainsi définie au cours du temps. Est-ce satisfaisant ? \n",
    "\n",
    "Refaire le calcul et les graphiques avec $N=200$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c54e2a69-acd3-41d7-8086-3940d2f25a8d",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "c0f2275a-17f0-4840-afb6-00a64e1dd5e3",
   "metadata": {},
   "source": [
    "C'est fini pour les TPs !"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9c30b731-0740-430d-b537-fa92fd63e292",
   "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.11.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}