{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Quelques méthodes d'intégration numérique"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Marc Lorenzi - 15 avril 2018"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from math import *\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Le but de ce notebook est de comparer quelques méthodes d'intégration numérique. Nous nous intéresserons à trois métodes : les méthodes des rectangles et des trapèzes vues en cours, et la méthode de Gauss. Pas de théorie ici, juste des illustrations de la théorie ou des conjectures pour la méthode de Gauss."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1. La méthode des rectangles"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Soit $f:[a,b]\\to\\mathbb R$. Pour tout $n\\ge 1$, soit $S_n=\\frac{b-a}{n}\\sum_{k=0}^{n-1}f(a + k\\frac{b-a}{n})$. Si $f$ est continue alors $S_n\\to\\int_a^b f$ lorsque $n$ tend vers l'infini. La raison essentielle à cela est que toute fonction continue peut être approchée de façon \"convenable\" par des fonctions en escalier. L'aire sous la courbe peut donc être approchée par des sommes d'aires de rectangles."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def rectangles(f, a, b, n):\n",
    "    d = (b - a) / n\n",
    "    s = 0\n",
    "    for k in range(n):\n",
    "        s += f(a + k * d)\n",
    "    return s * d"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rectangles(exp, 0, 1, 1000)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Si $f$ est lipschitzienne, la théorie nous dit que l'erreur commise en approchant l'intégrale de $f$ par $S_n$ est un $O(\\frac 1 n)$. Un petit tracé de l'erreur avec des axes logarithmiques permet d'illustrer cela. En abscisse, le nombre de points, $n$. En ordonnée, l'erreur commise. On choisit d'intégrer l'exponentielle entre 0 et 1, mais rien ne vous empêche de tenter autre chose. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def erreur(methode, f, a, b, n, vraie_valeur):\n",
    "    return abs(methode(f, a, b, n) - vraie_valeur)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "erreur(rectangles, exp, 0, 1, 1000, exp(1) - 1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ks = [10 * k for k in range(1, 1001)]\n",
    "srect = [erreur(rectangles, exp, 0, 1, k, exp(1) - 1) for k in ks]\n",
    "plt.loglog(ks, srect)\n",
    "plt.grid()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Retrouve-t-on bien cette fameuse erreur en $O(\\frac 1 n)$ ? Il semble que oui mais il y a mieux à tracer. Appelons $E_n$ l'erreur commise lorsqu'on intègre avec $n$ points. Imaginons que $E_n\\sim \\frac c {n^\\alpha}$ où $c>0$ et $\\alpha$ est le nombre qui nous intéresse. Nous soupçonnons que $\\alpha = 1$. On a $\\frac{E_{n}}{E_{2n}}\\sim 2^\\alpha$ et donc tend vers $2^\\alpha$ lorsque $n$ tend vers l'infini. Traçons donc $\\lg\\frac{E_{n}}{E_{2n}}$ où $\\lg$ désigne le logarithme en base 2."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def estimation_erreur(methode, f, a, b, n, vraie_valeur):\n",
    "    e1 = erreur(methode, f, a, b, n, vraie_valeur)\n",
    "    e2 = erreur(methode, f, a, b, 2 * n, vraie_valeur)\n",
    "    return log(e1 / e2) / log(2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "estimation_erreur(rectangles, exp, 0, 1, 1000, exp(1) - 1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ks = [10 * k for k in range(1, 101)]\n",
    "srect = [estimation_erreur(rectangles, exp, 0, 1, k, exp(1) - 1) for k in ks]\n",
    "plt.plot(ks, srect)\n",
    "plt.grid()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "On retrouve bien $\\alpha = 1$. Existe-t-il des fonctions pour lesquelles l'erreur commise est meilleure que prévu ? Eh bien oui. Prenons la fonction sinus sur $[0,\\pi]$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ks = [10 * k for k in range(1, 101)]\n",
    "srect = [estimation_erreur(rectangles, sin, 0, pi, k, 2) for k in ks]\n",
    "plt.plot(ks, srect)\n",
    "plt.grid()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "On voit que l'erreur est un $O(\\frac 1 {n^2})$. Il y a des raisons subtiles à cela, que nous ne détaillerons pas ici. Mais existe-t-il une méthode qui donne à tous les coups (pour des fonctions suffisamment régulières) une erreur en $O(\\frac 1 {n^2})$ ? Oui, par exemple la méthode des trapèzes."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2. La méthode des trapèzes"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Toute fonction suffisamment régulière peut être approchée de façon \"très convenable\" par des fonctions affines par morceaux. Il s'ensuit que l'aire sous la courbe peut être approchée par des sommes d'aires de trapèzes. Le cours montre que ces sommes sont à peu près celles de la méthode des rectangles, à un terme près qui est $\\frac 1 {2n}(b-a)(f(b)-f(a))$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def trapezes(f, a, b, n):\n",
    "    r = rectangles(f, a, b, n)\n",
    "    return r + (b - a) * (f(b) - f(a)) / (2 * n)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "trapezes(exp, 0, 1, 1000)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ks = [10 * k for k in range(1, 101)]\n",
    "srect = [estimation_erreur(trapezes, exp, 0, 1, k, exp(1) - 1) for k in ks]\n",
    "plt.plot(ks, srect)\n",
    "plt.grid()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Pour une fonction suffisamment régulière ($\\mathcal C^2$ par exemple), l'erreur commise est cette fois-ci un $O(\\frac 1 {n^2})$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3. Méthode de Gauss"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Peut-on faire mieux que les trapèzes ? Oui. Nous nous intéressons ici à une méthode due à Gauss. Pour être exacts, il s'agit de toute une famille de méthodes dont nous n'examinerons qu'un cas particulier. L'idée est d'approximer l'intégrale $\\int_a^b f$ par des sommes $\\sum_{k=0}^{n-1} a_k f(x_k)$ où les $x_k$ sont choisis très judicieusement dans le segment $[a,b]$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3.1 Un cas particulier"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Soit $f:[a,b]\\to\\mathbb R$ une fonction \"régulière\". Posons $m=\\frac{a+b}{2}$, $d=\\frac{b-a}{2}$ et $r=\\sqrt{\\frac 3 5}$. Considérons $S=\\frac 5 9 f(m-r d)+\\frac 8 9f(m)+\\frac 5 9 f(m+r d)$. On peut montrer que $S$ est une bonne approximation (en un sens que nous ne préciserons pas ici) de $\\int_a^b f$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def gauss0(f, a, b):\n",
    "    r = sqrt(3. / 5.)\n",
    "    m = (a + b) / 2.\n",
    "    d = (b - a) / 2.\n",
    "    return d * (5 * f(m - d * r) + 8 * f(m) + 5 * f(m + d * r)) / 9"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gauss0(exp, 0, 1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gauss0(exp, 0, 1) - exp(1) + 1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Une erreur de $10^{-6}$, alors qu'on a additionné 3 termes. C'est prometteur."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3.2 Cas général"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "La formule générale d'intégration de Gauss consiste à découper le segment $[a,b]$ en $n$ segments de tailles égales, et à appliquer le cas particulier sur chacun des segments."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def gauss(f, a, b, n):\n",
    "    s = 0\n",
    "    d = (b - a) / n\n",
    "    for k in range(n):\n",
    "        s = s + gauss0(f, a + k * d, a + (k + 1) * d)\n",
    "    return s"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "L'approximation obtenue est excellente."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "erreur(gauss, exp, 0, 1, 100, exp(1) - 1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "On obtient une erreur de $10^{-16}$ alors qu'avec la méthode des rectangles l'erreur était de $\\frac 1 {100}$ et avec la métode des trapèzes de $10^{-4}$. Comme précédemment, traçons l'erreur commise en fonction de $n$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ks = [10 * k for k in range(1, 1001)]\n",
    "sgauss = [erreur(gauss, exp, 0, 1, k, exp(1) - 1) for k in ks]\n",
    "plt.loglog(ks,sgauss)\n",
    "plt.grid()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Ah mais c'est quoi ça ? Eh bien voilà ce qui arrive quand on utilise des outils trop puissants ! En fait, l'erreur commise par la méthode de Gauss devient rapidement inférieure aux erreurs d'arrondis. Ce que nous traçons pour les \"grandes\" valeurs de $n$ ce n'est pas l'erreur commise, c'est du bruit. Au vu de la courbe ci-dessus, pour l'exemple de l'exponentielle entre 0 et 1, il est inutile de prendre une valeur de $n$ supérieure à ... combien au juste ? Retraçons le graphe avec une échelle semi-logarithmique et des valeurs de $k$ inférieures à 100."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ks = [k for k in range(1, 101)]\n",
    "sgauss = [erreur(gauss, exp, 0, 1, k, exp(1) - 1) for k in ks]\n",
    "plt.semilogy(ks,sgauss)\n",
    "plt.grid()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Inutile de dépasser 40 points pour intégrer. Mais quel est l'ordre de grandeur de l'erreur commise par la méthode de Gauss ? "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ks = [k for k in range(1, 15)]\n",
    "srect = [estimation_erreur(gauss, exp, 0, 1, k, exp(1) - 1) for k in ks]\n",
    "plt.plot(ks, srect)\n",
    "plt.grid()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$\\alpha = 6$ semble prometteur. Tout ceci reste bien entendu à prouver, et il s'avère que c'est bien la bonne valeur. Pour terminer, il faut noter que la régularité de la fonction à intégrer est importante pour que l'erreur soit petite. Essayons avec la fonction racine carrée sur 0, 1. On a $\\int_0^1\\sqrt x\\,dx=\\frac 2 3$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ks = [k for k in range(1, 15)]\n",
    "srect = [estimation_erreur(gauss, sqrt, 0, 1, k, 2./3.) for k in ks]\n",
    "plt.plot(ks, srect)\n",
    "plt.grid()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "On voit qu'ici l'erreur commise est un $O(1 / n^{\\frac 3 2})$, bien loin du $O(\\frac{1}{n^{6}})$ espéré."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "J'espère que vous avez compris le principe : ex-pé-ri-men-tez, c'est le grand intérêt des notebooks interactifs. Allez, une dernière question : quel est l'ordre de grandeur de l'erreur commise lorsqu'on intègre $\\sqrt x$ avec la méthode des rectangles ? Des trapèzes ? Lorsqu'on intègre $\\sin x$ avec la méthode de Gauss ? À vous de poursuivre ..."
   ]
  },
  {
   "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.6.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}