{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Exercice: Formule du rectangle et du trapèze" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "import matplotlib.pyplot as plt\n", "import numpy as np" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "On considère le calcul de l'intégrale\n", "$$\n", "I=\\int_{0}^1 f(x) \\, dx,\n", "$$\n", "où $f(x)$ est une fonction continue sur $[0,1]$.\n", "\n", "1. Ecrivez une fonction `midpoint` qui implémente la formule composite du rectangle (point\n", "milieu) pour l'approximation de l'intégrale ci-dessus. Pour permettre le choix d'un intervalle d'intégration de la forme $[a,b]$, le nombre de sous-intervalles $N$\n", "et la fonction $f(x)$, définie par la commande `f = lambda x : ... `\n", " utilisez la structure suivante:\n", "\n", "```python\n", "def midpoint( a, b, N, f ) :\n", " # [a,b] : inteval\n", " # N : number of subintervals\n", " # f : fonction to integrate using the midpoint rule\n", "``` \n", " \n", "2. Testez le code en considérant la fonction $f(x)=x^2$ et\n", "$M=10$. Tracez ensuite le graphe de l'erreur $|I(f)-I_{\\mathrm{\\textit{pm}}}^c (f)|$\n", "en fonction de $N$ (utilisez à cette fin les commandes `plt.xscale('log'); plt.yscale('log')`), en choisissant $a=0$, $b=1$,\n", "et $M = 10^1, 10^2, 10^3, \\ldots 10^5$. Quel est l'ordre de\n", "convergence de la formule composite du rectangle par rapport à $H=(b-a)/M$? Donnez une\n", "interprétation des résultats d'après la théorie.\n", "\n", "3. Modifiez le code du point 1. pour permettre le calcul\n", "de l'intégrale avec la formule composite du trapèze.\n", "Tracez le graphe de l'erreur $|I(f)-I_{\\mathrm{\\textit{t}}}^c (f)|$\n", "pour les mêmes valeurs de $M$.\n", "Comparez les résultats avec ceux obtenus\n", "au point 2.\n", "\n", "```python\n", "def trapeziodal( a, b, N, f ) :\n", " # [a,b] : inteval\n", " # N : number of subintervals\n", " # f : fonction to integrate using the trapezoidal rule\n", "``` \n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Partie 1" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def midpoint( a, b, N, f ) :\n", " # [a,b] : inteval\n", " # M : number of subintervals\n", " # f : fonction to integrate using the midpoint rule\n", "\n", " # (completez)\n", "\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Partie 2" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "f = lambda x : x**2\n", "\n", "a = 0; b = 1; N = 10;\n", "\n", "intmp = midpoint( a, b, N, f )\n", "\n", "print(intmp)\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "En sachant que la valeur exacte de l'intégrale est $I(f)=1/3$,\n", "on peut écrire les commandes suivantes:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "N = 10**np.linspace(1,5,5).astype(int)\n", "\n", "errmp = []\n", "for i in range(5) :\n", " # (completez)\n", "\n", "# (completez)\n", " \n", "\n", " \n", " \n", "slope = ( np.log(errmp[4]) - np.log(errmp[0]) ) / ( np.log(N[4]) - np.log(N[0]) ) \n", "\n", "print(f'La pente est de {slope:.6f}')\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Donc, l'erreur décroît comme la puissance -2 du paramètre $N$.\n", "Par conséquent, l'ordre de convergence par rapport à $H = 1/N$ est $2$.\n", "\n", "\n", "#### Partie 3 \n", "\n", "On peut modifier le code du point milieu pour\n", "implémenter la méthode du trapèze comme suit." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def trapezoidal( a, b, N, f ) :\n", " # [a,b] : inteval\n", " # M : number of subintervals\n", " # f : fonction to integrate using the trapezoidal rule\n", "\n", " # (completez)\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "errtrap = []\n", "\n", "# (completez)\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "On observe que la précision est encore $\\mathcal{O}(h^2)$, mais\n", "que l'erreur est plus grande (le double, précisément)\n", "que dans le cas du point milieu.\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Exercice: Formule de Simpson" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "import matplotlib.pyplot as plt\n", "import numpy as np" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "On considère la fonction $f \\, : \\, [a,b] \\rightarrow \\mathbb R$ dans $C^0([a,b])$;\n", "on est intéressé à approcher l'intégrale $I(f) = \\int_a^b f(x) \\, dx$. \n", "\n", "1. Ecrivez une fonction `Simpson` qui implémente la formule composite de Simpson pour l'approximation de l'intégrale ci-dessus. \n", "Utilisez la structure suivante:\n", "\n", "```python\n", "def Simpson( a, b, N, f ) :\n", " # [a,b] : inteval\n", " # N : number of subintervals\n", " # f : fonction to integrate using the Simpson rule\n", "```\n", "\n", "2. Testez ensuite pour quels monômes $f(x) = x^d$ cette formule intègre exactement la fonction, pour $d=0,1,...$ sur l'intervalle $[1,4]$, avec $N=1$ et ensuite $N=10$.\n", " \n", "3. Vérifiez numériquement pour quelques polynômes que la fonction ainsi écrite est linéaire en $f$ pour $N=10$\n", "\n", "*Suggestion*: En utilisant une fonction `lambda` il est possible de décider le paramètre $d$ d'un monôme à un moment ultérieur:\n", "\n", "```python\n", "monomial = lambda x : x**d\n", "```\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Partie 1" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# This function just provides the Simpson quadrature rule.\n", "def Simpson( a, b, N, f ) :\n", " # [a,b] : inteval\n", " # N : number of subintervals\n", " # f : fonction to integrate using the trapezoidal rule\n", "\n", " # (completez)\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Partie 2\n", "\n", "$$\\int_1^4 1 dx= 3$$\n", "$$\\int_1^4 x^d dx = \\frac1{d+1} x^{d+1} \\vert_1^4 = \\frac{4^{d+1} - 1}{d+1}$$" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Checking simpson fonction\n", "a = 1; b = 4;\n", "\n", "# with lambda fonctions, it is possible to determine a parameter (here d)\n", "# at a later moment\n", "monomial = lambda x : x**d\n", "\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "N = 1\n", "for d in range(7) :\n", " # (completez)\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "N = 10\n", "\n", "for d in range(7) :\n", " # (completez)\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Partie 3\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# (completez)\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Exercice: Formule du rectangle, du trapèze et de Simpson" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "On considère la fonction $f \\, : \\, [a,b] \\rightarrow \\mathbb R$ dans $C^0([a,b])$;\n", "on est intéressé à approcher l'intégrale $I(f) = \\int_a^b f(x) \\, dx$. \n", "\n", "Plus précisément on prend\n", "$f(x) = \\sin( \\frac72 \\, x ) + e^x- 1$ avec $a=0$ et $b=1$ ($f \\in C^{\\infty}([a,b])$) et on\n", "peut calculer $I(f) = \\frac27 \\, (1-\\cos(\\frac72 )) + e - 2$.\n", "\n", "1. Calculez une approximation de l'intégrale en utilisant les formules du rectangle, du trapèze et\n", "de Simpson *simples*, c-à-d avec un seul intervalle. \n", "\n", "2. Calculez une approximation de l'intégrale en utilisant les fonctions `Midpoint`, `Trapezoidal` et `Simpson` (déjà codées).\n", " avec $N=10$ sous-intervalles de même taille.\n", " On notera les valeurs approchées de l'intégrale\n", " $I_{mp}^c(f)$, $I_{t}^c(f)$, and $I_{s}^c(f)$, respectivement.\n", "\n", "3. Répétez le point 2. avec $N=2^{k}$ pour $k=2,\\ldots,7$ et calculez les erreurs $E^c_{mp}(f) := | I(f) - I_{mp}^c(f)|$, $E^c_{t}(f) := | I(f) - I_{t}^c(f)|$, et $E^c_{s}(f) := | I(f) - I_{s}^c(f)|$.\n", " Dessinez les erreurs en fonction de $H=(b-a)/N$ sur une échelle logarithmique sur les deux axes.\n", " Quel est l'ordre de convergence de ces méthodes ? Est-ce en accord avec la théorie ?\n", " Motivez votre réponse.\n", "\n", "4. On prend maintenant $f(x) = x^d$, $a=0$ et $b=1$, avec $d \\in \\mathbb N$. L'intégrale de $f$ vaut $I(f) = 1/(d+1)$.\n", " Vérifiez numériquement les degrés d'exactitude de chacune des formules de quadrature du point 1.\n", " Pour cela, il faut choisir plusieurs valeurs de $d = 0,1,2,\\ldots$. \n", " Motivez votre réponse.\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "import matplotlib.pyplot as plt\n", "import numpy as np\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Partie 1" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "f = lambda x : np.sin(7/2*x) + np.exp(x) - 1\n", "\n", "a = 0; b = 1\n", "Iexact = 2/7*(1-np.cos(7/2) ) + np.exp(1) - 2\n", "\n", "x = np.linspace(a,b,1000)\n", "y=f(x)\n", "\n", "plt.plot(x, y, 'b')\n", "plt.show()" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "intMD = \n", "intTrap = \n", "intSimp = \n", "\n", "print(f'exact \\t rect \\t\\t trap \\t\\t Simpson')\n", "print(f'{Iexact:.4f} \\t {intMD:.4f} \\t {intTrap:.4f} \\t {intSimp:.4f}')\n", " " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Partie 2" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "N = 10\n", "intMD = Midpoint(a,b,N,f)\n", "intTrap = Trapezoidal(a,b,N,f)\n", "intSimp = Simpson(a,b,N,f)\n", "\n", "print(f'exact \\t rect \\t\\t trap \\t\\t Simpson')\n", "print(f'{Iexact:.4f} \\t {intMD:.4f} \\t {intTrap:.4f} \\t {intSimp:.4f}')\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Partie 3" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "errmp = []\n", "errtrap = []\n", "errSim = []\n", "\n", "# (completez)\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Partie 4\n", "\n", " The function $f(x)=x^d$ is a polynomial of degree $d$ for $d \\in \\mathbb N$. The simple midpoint, trapezoidal, and Simpson quadrature formulas possesses degree of exactness equal to $1$, $1$, and $3$, respectively." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "N = 1\n", "# with lambda fonctions, it is possible to determine a parameter (here d)\n", "# at a later moment\n", "monomial = lambda x : x**d\n", "\n", "# (completez)\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# (completez)\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# (completez)\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "De ces simulations nous vérifions que les monômes $f(x)=x^d$ sont intégrées exactement pour les degrées $d=0$ et $d=1$ dans le cas des formules du rectangle et du trapèze, et pour $d=0,1,2,3$ dans le cas de la formule de Simpson." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Exercice: Formule du rectangle et du trapèze" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "import matplotlib.pyplot as plt\n", "import numpy as np\n", "\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "On considère, sur l'intervalle $[-1,1]$, la fonction $f$ suivante:\n", "\n", "$$f(x) = \\left\\{ \\begin{array}{cc}\n", "\\frac{}{}e^x & \\mbox{ si } x\\leq 0 \\\\\n", "1 & \\mbox{ si } x > 0 \\\\\n", "\\end{array} \\right.$$\n", "\n", "On peut définir une telle fonction en utilisant la commande\n", "\n", "```python\n", "f = lambda x : np.exp(x)*(x<=0) + (1)*(x>0)\n", "```\n", "\n", "\n", "1. Utilisez $1000$ points équirépartis dans l'intervalle $[-1,1]$ pour afficher la fonction $f$ (utilisez la commande `axis` pour recadrer l'image).\n", "\n", "2. On s'intéresse à présent à l'intégrale $I=\\int_{-1}^1 f(x) \\ dx$. On peut calculer la valeur de $I$ analytiquement et on trouve $I=2-\\frac{1}{e} \\cong 1.6321$. Calculez des valeurs approchées de $I$ en considérant les formules du point milieu, du trapèze et de Simpson avec $N=1,9,99,999$ où $N$ est le nombre de sous-intervalles de formules composites. Utilisez les fonctions `midpoint`, `trapezoidal` et `simpson` (déjà codées).\n", "\n", "3. Reportez les erreurs calculées au point (b) dans un graphe montrant l'erreur en fonction de $H$ avec des échelles logarithmiques.\n", "\n", "4. Estimez, à partir des graphes obtenus au point précédent, l'ordre de chacune des méthodes. Comparez-les avec les ordres donnés au cours.\n", "Y a-t-il des différences? Pourquoi? (Regardez les dérivées de $f$ ).\n", "\n", "5. Pourquoi obtient-on de bien meilleurs résultats pour la méthode de Simpson avec un nombre pair d'intervalles qu'avec un nombre impair (essayez avec $N=99$ et ensuite $N=100$ )?\n", "\n", "6. Refaites l'exercice avec $N=2,10,100,1000$.\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Partie 1" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "f = lambda x : np.exp(x)*(x<=0) + (1)*(x>0)\n", "\n", "a = -1; b = 1\n", "x = np.linspace(a,b,1000)\n", "y=f(x)\n", "\n", "plt.plot(x, y, 'b')\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Partie 2\n", "\n", "On calcule tout d'abord la valeur exacte de l'intégrale, puis on calcule les erreurs, que l'on stocke dans des vecteurs:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "Iexact = 2 - np.exp(-1)\n", "\n", "# (completez)\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Partie 3" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# (completez)\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Les graphes pour toutes les méthodes sont des droites. On remarque que lorsque $H$ est divisé par $10$, les \n", "erreurs sont environ divisées par $100$. On peut donc supposer que les ordres sont $2$ par rapport à $H$. Pour le confirmer,\n", "on peut ajouter sur le graphe la pente représentant l'ordre $2$ et vérifier qu'elle est parallèle aux\n", "droites d'erreur." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Partie 4" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# (completez)\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Les graphiques sont parallèles à la droite $H^2$ et donc l'ordre est $2$ pour toutes les méthodes. Pour la\n", "méthode du point milieu et du trapèze, c'est l'ordre auquel on s'attend d'après la théorie. Par contre, pour Simpson, on\n", "pouvait s'attendre à un ordre $4$. Le problème vient du fait que la fonction n'est pas très régulière, puisque \n", "$f \\in C^0([-1,1])$ mais $f \\not\\in C^1([-1,1])$ puisque sa dérivée n'est pas continue en $x=0$. Dans la théorie, on demande\n", "$f \\in C^4([-1,1])$ pour assurer l'ordre $4$. On obtient tout de même l'ordre $2$ car, mis à part en $x=0$, la fonction $f$\n", "est très régulière." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Partie 5\n", "\n", "Si on regarde l'erreur pour la méthode de Simpson avec $N=99$ et $N=100$ sous-intervalles" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "print (np.abs( Simpson( a, b, 99, f) - Iexact ) )\n", "print (np.abs( Simpson( a, b, 100, f) - Iexact ) )\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "On obtient des valeurs d'environ $1.70 \\ 10^{-5}$ et $3.51 \\ 10^{-11}$. Les erreurs sont très différentes !\n", "Voilà pourquoi:\n", "\n", "* si $N$ est pair, ...\n", "\n", "* si $N$ est impair, ..." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Partie 6\n", "\n", "Il faut relancer les parties 2 et 3 avec `N= [2,10,100,1000]` à la place de `N = [1,9,99,999]`\n" ] }, { "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.7.6" } }, "nbformat": 4, "nbformat_minor": 4 }