{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Résolution numérique des équations différentielles"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Marc Lorenzi\n",
    "\n",
    "15 décembre 2020"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import math"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.rcParams['figure.figsize'] = (10, 10)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Dans tout ce notebook les fonctions seront définies sur un segment $[a,b]$ de $\\mathbb R$, à valeurs dans $\\mathbb R$ ou parfois dans $\\mathbb R^p$ ($p\\ge 2$).\n",
    "\n",
    "**Notations** : \n",
    "\n",
    "- On note la « variable » par la lettre $t$. \n",
    "- Lorsque la fonction est à valeurs réelles, on la note par une lettre minuscule du genre $x, u, v$, etc. Par exemple, si l'on résout l'équation différentielle $x'=x+t$, on notera une solution de l'équation par la lettre $x$, l'image du réel $t$ par la fonction $x$ est notée $x(t)$, etc.\n",
    "- Les fonctions à valeurs dans $\\mathbb R^p$ sont des fonctions à valeurs vectorielles, on les notera de préférence par des lettres majuscules, $X, Z, etc$.\n",
    "\n",
    "Les fonctions seront aussi régulières que nécessaire."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1. Introduction"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1.1 Équation d'ordre 1 : l'idée générale"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Soit \n",
    "\n",
    "$$(\\mathcal E)\\quad x'=F(t,x)$$\n",
    "\n",
    "une équation différentielle linéaire d'ordre 1 « résolue en $x'$ ». Soit $x_0\\in\\mathbb R$. Soient $a,b\\in\\mathbb R$, $a< b$. Sous certaines conditions de régularité sur la fonction $F$ que je n'expliciterai pas ici, il existe une unique solution $x:[a,b]\\to\\mathbb R$ de l'équation différentielle sur le segment $[a,b]$ telle que $x(a)=x_0$.\n",
    "\n",
    "Supposons connue une **méthode** qui nous permet, pour $t\\in[a,b]$ et $h$ « petit », de calculer à partir d'une valeur approchée de $x(t)$, une valeur approchée de $x(t+h)$. Si l'on connaît $x(a)$, une simple boucle permet alors de calculer une valeur approchée de $x$ en les points d'une subdivision de $[a,b]$ de pas $h$, en appliquant itérativement cette méthode.\n",
    "\n",
    "Tout réside dans le sens du mot « approchée ». Quelle est la précision de __l'approximation__ ? Tout dépend évidemment de la __méthode__. Nous pouvons déjà écrire une fonction générale de résolution approchée d'une équation différentielle d'ordre 1 ... "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "La fonction `dsolve` prend en paramètres :\n",
    "\n",
    "- Une fonction $F$, qui est le « $F$ » de l'équation différentielle à résoudre.\n",
    "- Deux réels $a$ et $b$ qui sont les bornes de l'intervalle sur lequel on résout.\n",
    "- Un réel $x_0$, la valeur initiale de notre solution en $a$.\n",
    "- Une entier $n$ qui est le nombre de points de $[a,b]$ en lesquels nous voulons une valeur approchée de la solution.\n",
    "- Une fonction `method`, qui est ... la __méthode__ permettant de calculer $x(t+h)$ à partir de $x(t)$, connaissant $x(t)$, $F$, $t$ et $h$. Un appel à  `method(F, t, x, h)` renvoie une valeur approchée de $x(t+h)$ lorsque le paramètre `x` est une valeur approchée de $x(t)$.\n",
    "\n",
    "La fonction `dsolve` renvoie un couple $(ts, xs)$ où $ts$ est la liste des points d'une subdivision régulière de $[a,b]$ de pas $\\frac{b-a}n$ et $xs$ est la liste des valeurs approchées de la solution $x$ aux points de la subdivision. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "__Remarque__ : les professionnels ne disent pas « méthode » mais __schéma__. Étant un amateur, nous continuerons à employer le mot « méthode »."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def subdivision(a, b, n):\n",
    "    h = (b - a) / n\n",
    "    return [a + k * h for k in range(n + 1)]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def dsolve(F, a, b, x0, n, method):\n",
    "    h = (b - a) / n\n",
    "    ts = subdivision(a, b, n)\n",
    "    xs = [x0]\n",
    "    x = x0\n",
    "    for k in range(n):\n",
    "        t = ts[k]\n",
    "        x = method(F, t, x, h)\n",
    "        xs.append(x)\n",
    "    return (ts, xs)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Testons .. avec une méthode absurde :-), juste pour voir si la fonction fonctionne. Résolvons $x'=x$ sur $[0, 1]$ avec $x(0)=1$ et une subdivision de 10 points. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "F = lambda t, x: -x\n",
    "method = lambda F, t, x, h: 0\n",
    "print(dsolve(F, 0, 1, 1, 10, method))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "La solution cherchée est l'exponentielle. Évidemment la liste $[1,0,0,\\ldots,0]$ est une approximation assez effroyable des valeurs de $x\\mapsto e^x$ aux points $0,\\frac{1}{10},\\frac{2}{10},\\ldots, 1$. Mais on s'en doutait. On a évidemment compris que tout réside dans la __méthode__. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1.2 Équation d'ordre 2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Nous allons bientôt voir pourquoi il est utile de considérer des équations différentielles \n",
    "\n",
    "$$X'=G(t, X)$$\n",
    "\n",
    "où $X:[a,b]\\to\\mathbb R^p$ est une fonction à valeurs __vectorielles__ dans $\\mathbb R^p$. Une telle fonction $X$ est définie par une égalité du type \n",
    "\n",
    "$$X(t)=(x_1(t),\\ldots,x_p(t))$$\n",
    "\n",
    "où les $x_i:[a,b]\\to\\mathbb R$ sont des fonctions à valeurs réelles. Dériver $X$, c'est tout simplement dériver chacune des fonctions $x_i$ : \n",
    "\n",
    "$$X'(t)=(x'_1(t),\\ldots,x'_p(t))$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "S'intéresser à des fonctions à valeurs vectorielles c'est juste pour faire _tendance_ ? Non, pas du tout ! Considérons une équation différentielle d'ordre 2, \n",
    "\n",
    "$$x''=F(t,x,x')$$\n",
    "\n",
    "Posons $X=(x,x')$. Posons également \n",
    "\n",
    "$$G(t, (u, v)) = (v, F(t, u, v))$$\n",
    "\n",
    "Alors \n",
    "\n",
    "$$\\begin{array}{lll}\n",
    "X'&=&(x', x'')\\\\\n",
    "&=&(x', F(t,x,x'))\\\\\n",
    "&=&G(t,(x, x'))\\\\\n",
    "&=&G(t,X)\n",
    "\\end{array}$$\n",
    "\n",
    "Bref, résoudre l'équation différentielle d'ordre 2, $x''=F(t,x,x')$ revient à résoudre l'équation vectorielle d'ordre 1 : \n",
    "\n",
    "$$X'=G(t,X)$$\n",
    "\n",
    "Bien sûr, ceci peut se généraliser à des équations différentielles d'ordre 3, 4, etc. Par exemple, pour une équation d'ordre 3 on prendrait $X=(x,x',x'')$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "La fonction `dsolve2` ci-dessous permet de \"résoudre\" des équations d'ordre 2. Ici, $x_0$ et $x_1$ sont respectivement les valeurs de la fonction $x$ et de sa dérivée $x'$ au point $a$. Les autres paramètres ont la même signification que pour la fonction `dsolve`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def dsolve2(F, a, b, x0, x1, n, method):\n",
    "    X0 = [x0, x1]\n",
    "    G = lambda t, X: (X[1], F(t, X[0], X[1]))\n",
    "    ts, Xs = dsolve(G, a, b, X0, n, method)\n",
    "    return (ts, [X[0] for X in Xs])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Le seul souci que pour l'instant on ne peut résoudre aucune équation vu qu'il nous manque des méthodes. Mais c'est pour ça qu'on est en prépa, pour acquérir des méthodes, non ? Nous allons en voir deux :\n",
    "\n",
    "- La __méthode d'Euler__ est très facile à implémenter, intuitivement claire, ... et très inefficace.\n",
    "- La __méthode de Runge-Kutta__ est un peu plus délicate à écrire, pas du tout intuitive ... et terriblement efficace."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2. Zone de repos ..."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Nous aurons besoin sous peu d'additionner des vecteurs de $\\mathbb R^p$ et de les multiplier par un réel. Nous représenterons ces vecteurs en Python par des listes de longueur $p$. Écrivons des fonctions d'addition et de multiplication par un réel. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def mul(X, mu):\n",
    "    return [x * mu for x in X]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mul([1, 3, 2], 5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "__Question__ : pourquoi `mu` et pas `lambda` ???"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def add(X1, X2):\n",
    "    return [X1[k] + X2[k] for k in range(len(X1)) ]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "add([1, 2, 6], [3, 2, 4])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3. Euler"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3.1 La méthode, enfin"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Considérons l'équation différentielle \n",
    "\n",
    "$$x'=F(t,x)$$\n",
    "\n",
    "Soit $x$ une fonction solution de cette équation sur un intervalle $I$. Soit $t\\in I$. Pour $h$ réel tel que $x+h\\in I$, on a"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$x(t+h) = x(t)+\\int_t^{t+h}x'(u)\\,du=x(t)+\\int_t^{t+h}F(u,x(u))\\,du$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Si $h$ est « petit », une approximation très crue de l'intégrale est \n",
    "\n",
    "$$\\int_t^{t+h}F(u,x(u))\\,du\\simeq hF(t,x(t))$$\n",
    "\n",
    "On peut donc espérer (il faut être optimistes) que \n",
    "\n",
    "$$x(t+h)\\simeq x(t)+hF(t,x(t))$$\n",
    "\n",
    "Je n'irai pas plus avant dans la théorie, ce n'est pas le but de ce notebook. Mais nous tenons là une méthode. Est-elle efficace ?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def euler_method(F, t, x, h):\n",
    "    dx = F(t, x) * h\n",
    "    return x + dx"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3.2 Itération"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def euler(F, a, b, x0, n):\n",
    "    return dsolve(F, a, b, x0, n, euler_method)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Résolvons l'équation $y'=y$ sur $[0,2]$ avec une valeur initiale de 1 en 0. La solution exacte est bien sûr la fonction exponentielle."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ts, xs = euler(lambda t, x:x, 0, 2, 1, 100)\n",
    "print(xs)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Traçons sur un même graphe la solution approchée et la solution exacte."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ts, xs = euler(lambda t, x:x, 0, 2, 1, 100)\n",
    "plt.plot(ts, xs, 'k')\n",
    "plt.plot(ts, [math.exp(t) for t in ts], 'r')\n",
    "plt.grid()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Quelle est l'erreur commise ?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "err = [abs(xs[k] - math.exp(ts[k])) for k in range(len(ts))]\n",
    "plt.plot(ts, err, 'k')\n",
    "plt.grid()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Une erreur de $0.14$ à la dernière étape, ce n'est pas très bon, c'est le moins qu'on puisse dire. On peut montrer que la méthode d'Euler est ce que l'on appelle une __méthode d'ordre 1__ : si on appelle $E_k$ l'erreur commise à l'étape $k$, on a pour tout $k$, **SI** $h$ est assez petit, \n",
    "\n",
    "$$E_k\\le Ch$$\n",
    "\n",
    "où $C$ est une « constante » et $h$ est le pas de la méthode (le dernier paramètre de `euler_method`). Cette « constante » $C$ dépend en fait de tout sauf de $h$. Elle dépend en particulier exponentiellement de la longueur de l'intervalle où l'on résout l'équation différentielle."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2.3 Fonctions à valeurs vectorielles"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "On recopie ici ce que l'on a fait pour les fonctions à valeurs réelles, sauf que l'on met des `mul` et des `add` à la place des multiplications et des additions. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def eulerV_method(G, t, X, h):\n",
    "    dX = mul(G(t, X), h)\n",
    "    return add(X, dX)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def eulerV(G, a, b, X0, n):\n",
    "    return dsolve(G, a, b, X0, n, eulerV_method)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Testons sur un petit exemple sans prétention. Considérons le __système différentiel__ \n",
    "\n",
    "$$(u', v')=(-v, u)$$\n",
    "\n",
    "où $u$ et $v$ sont deux fonctions inconnues de la variable $t$. Prenons les conditions initiales $u(0)=1$, $v(0)=0$. Que valent $u$ et $v$ ? Clairement, \n",
    "\n",
    "$$u(t)=\\cos t, \\ v(t)=\\sin t$$\n",
    "\n",
    "est une solution, et on peut montrer que c'est la seule. Si l'on trace la __courbe paramétrée__  $t\\mapsto (u(t),v(t))$ dans le plan, on obtient ... le cercle unité. Enfin si tout va bien."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def G(t, X): return [-X[1], X[0]]\n",
    "ts, Xs = eulerV(G, 0, 4 * math.pi, [1, 0], 1000)\n",
    "us = [X[0] for X in Xs]\n",
    "vs = [X[1] for X in Xs]\n",
    "plt.plot(us, vs, 'k')\n",
    "plt.grid()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "On voit ici que l'erreur en $O(h)$ est fatale. La périodicité des solutions fait éclater l'erreur au grand jour, notre cercle ne se referme pas. La méthode d'Euler est trop approximative pour permettre d'obtenir des valeurs approchées fiables des solutions d'une équation différentielle."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2.4 Équation différentielle d'ordre 2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "On a compris, Euler pas bien. Mais enfonçons le clou en résolvant l'équation $y''+y=0$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def euler2(F, a, b, x0, x1, n):\n",
    "    return dsolve2(F, a, b, x0, x1, n, eulerV_method)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ts, xs = euler2(lambda t, x, x1: -x, 0, 10 * math.pi, 0, 1, 1000)\n",
    "plt.plot(ts, xs, 'k')\n",
    "plt.plot(ts, [math.sin(t) for t in ts], 'r')\n",
    "plt.grid()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "En rouge, la vraie solution, qui est la fonction sinus. En noir, la solution approchée. Dessinons le graphe de l'erreur commise."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "err = [abs(xs[k] - math.sin(ts[k])) for k in range(len(ts))]\n",
    "plt.plot(ts, err, 'k')\n",
    "plt.grid()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "La méthode d'Euler permet d'obtenir rapidement des valeurs approchées de solutions d'équations différentielles, mais ces valeurs approchées sont de qualité plus que médiocre.\n",
    "\n",
    "Peut-on faire mieux ? Oui, et quitte à faire mieux, allons directement à __ce qui se fait de mieux__."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3. « La » méthode de Runge-Kutta"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3.1 Équation d'ordre 1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Il existe en fait __des__ méthodes de Runge-Kutta. Nous allons nous intéresser à la méthode dite d'ordre 4. Considérant toujours notre équation différentielle $x'=F(t, x)$, comment passer d'une valeur approchée de $x(t)$ à une valeur approchée de $x(t+h)$ ?\n",
    "\n",
    "On calcule tout d'abord successivement les 4 nombres :\n",
    "\n",
    "$$\\left\\{\\begin{array}{lll}\n",
    "k_0&=&h F(t,x)\\\\\n",
    "k_1&=&hF(t+\\frac h 2, x + \\frac {k_0} 2)\\\\\n",
    "k_2&=&hF(t+\\frac h 2, x + \\frac {k_1} 2)\\\\\n",
    "k_3&=& h F(t+h, x + k_2)\n",
    "\\end{array}\\right.$$\n",
    "\n",
    "On pose ensuite\n",
    "\n",
    "$$\\delta x = \\frac 1 6 (k_0 + 2k_1+2k_2+k_3)$$\n",
    "\n",
    "La méthode que nous examinons ici consiste à prendre \n",
    "\n",
    "$$x(t+h)\\simeq x + \\delta x$$\n",
    "\n",
    "Je ne donnerai ici aucune justification concernant l'excellence de cette méthode. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def rk4_method(F, t, x, h):\n",
    "    k0 = h * F(t, x)\n",
    "    k1 = h * F(t + 0.5 * h, x + 0.5 * k0)\n",
    "    k2 = h * F(t + 0.5 * h, x + 0.5 * k1)\n",
    "    k3 = h * F(t + h, x + k2)\n",
    "    dx = (1 / 6) * (k0 + 2 * k1 + 2 * k2 + k3)\n",
    "    return x + dx"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def rk4(F, a, b, x0, n):\n",
    "    return dsolve(F, a, b, x0, n, rk4_method)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Reprenons l'exemple de l'exponentielle."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ts, xs = rk4(lambda t, x:x, 0, 2, 1, 100)\n",
    "plt.plot(ts, xs, 'k')\n",
    "plt.plot(ts, [math.exp(t) for t in ts], 'r')\n",
    "plt.grid()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Tiens, on ne voit qu'une courbe ? Où est la courbe noire de la solution approchée ? Eh bien elle est __sous__ la courbe rouge, bien cachée ... tiendrait-on une bonne méthode ?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "err = [abs(xs[k] - math.exp(ts[k])) for k in range(len(ts))]\n",
    "plt.plot(ts, err, 'k')\n",
    "plt.grid()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Remarquez le $10^{-8}$ en haut à gauche du graphique : l'erreur maximale commise n'est pas 2, mais $2\\times 10^{-8}$. L'approximation obtenue est excellente."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3.2 Équation vectorielle"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Poursuivons nos tests. On réécrit la méthode pour des fonctions à valeurs vectorielles. Ce n'est qu'une réécriture avec des `add` et des `mul` de la fonction `rk4_method`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def rk4V_method(G, t, X, h):\n",
    "    K0 = mul(G(t, X), h)\n",
    "    K1 = mul(G(t + 0.5 * h, add(X, mul(K0, 0.5))), h)\n",
    "    K2 = mul(G(t + 0.5 * h, add(X, mul(K1, 0.5))), h)\n",
    "    K3 = mul(G(t + h, add(X, K2)), h)\n",
    "    dX = mul(add(add(K0, K3), mul(add(K1, K2), 2)), 1/6)\n",
    "    return add(X, dX)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def rk4V(G, a, b, X0, n):\n",
    "    return dsolve(G, a, b, X0, n, rk4V_method)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Retentons le tracé du cercle."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def G(t, X): return [-X[1], X[0]]\n",
    "ts, Xs = rk4V(G, 0, 4 * math.pi, [1, 0], 1000)\n",
    "us = [X[0] for X in Xs]\n",
    "vs = [X[1] for X in Xs]\n",
    "plt.plot(us, vs, 'k')\n",
    "plt.grid()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Cette fois-ci rien ne dépasse. Le cercle se referme."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3.3 Équation d'ordre 2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Rien à dire, copier-coller de ce qui a été fait pour la méthode d'Euler."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def rk42(F, a, b, x0, x1, n):\n",
    "    return dsolve2(F, a, b, x0, x1, n, rk4V_method)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Pour bien voir l'efficacité de la méthode, on intègre l'équation différentielle sur 20 périodes du sinus, avec 1000 pas."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ts, xs = rk42(lambda t, x, x1: -x, 0, 40 * math.pi, 0, 1, 1000)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.plot(ts, xs, 'k')\n",
    "plt.grid()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "err = [abs(xs[k] - math.sin(ts[k])) for k in range(len(ts))]\n",
    "plt.plot(ts, err, 'k')\n",
    "plt.grid()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "L'erreur maximale est donc environ $2\\times 10^{-4}$. Bien qu'il existe des méthodes encore meilleures, la méthode de Runge-Kutta d'ordre 4 est une méthode à envisager sérieusement lorsqu'on a à faire des simulations numériques."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "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.8.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}