diff --git a/3.2 Exercices.ipynb b/3.2 Exercices.ipynb index 5c91fff..e0d064b 100644 --- a/3.2 Exercices.ipynb +++ b/3.2 Exercices.ipynb @@ -1,730 +1,730 @@ { "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "## Intégration Numérique : Exercices" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# importing libraries used in this book\n", "import numpy as np\n", "import matplotlib.pyplot as plt" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Exercice Série 6, Ex 2 : Formule 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", "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", " M=3\n", " \n", " nodes = np.array([-1, 0, 1])\n", " weights = np.array([1./3., 4./3., 1./3.])\n", " \n", " # size of the subintervals\n", " ## H = << Complétez cette ligne\n", " # points defining intervals\n", " ## x = << Complétez cette ligne\n", " \n", " Lh = 0;\n", " \n", " z = np.zeros(M);\n", " \n", " for k in range(N) : \n", " # left of the subinterval, also first quadrature point\n", " z[0] = x[k]\n", " # right of the subinterval, also third quadrature point\n", " z[2] = x[k+1]\n", " # mid point, , also second quadrature point\n", " z[1] = (x[k] + x[k+1])/2\n", " # can also be computed as\n", " z = (x[k] + x[k+1])/2 + nodes*(x[k+1] - x[k])/2\n", "\n", " # local quadrature:\n", " ## Jgk = << Complétez cette ligne\n", "\n", " Lh = Lh + Jgk\n", "\n", " # approximate integral \n", " return H/2 * Lh" ] }, { "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", "# recording for which degrees the integral s exact (up to epsilon)\n", "exactDegree = -1\n", "epsilon = 1e-12" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "N = 1\n", "for d in range(7) :\n", " intsim = Simpson( a, b, N, monomial )\n", " intExact = (4**(d+1) - 1)/(d+1)\n", " print(f'Simpson on x^{d} : {intsim:.6f} - {intExact:.6f} = {intsim-intExact:.6e}')\n", " if np.abs(intsim-intExact) < epsilon :\n", " exactDegree = d\n", " \n", "print(f'Simpson with N = {N} is exact up to degree {exactDegree}')" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "N = 10\n", "exactDegree = -1\n", "\n", "for d in range(7) :\n", " intsim = Simpson( a, b, N, monomial )\n", " intExact = (4**(d+1) - 1)/(d+1)\n", " print(f'Simpson on x^{d} : {intsim:.6f} - {intExact:.6f} = {intsim-intExact:.6e}')\n", " if np.abs(intsim-intExact) < epsilon :\n", " exactDegree = d\n", " \n", "print(f'Simpson with N = {N} is exact up to degree {exactDegree}')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Partie 3\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": [ "maxD = 5;\n", "N = 10\n", "\n", "p = lambda x : np.polyval(coefs,x)\n", "\n", "# pre-computing integrls of monomials up to degree maxD\n", "intMono = np.zeros(maxD+1)\n", "for d in range(maxD+1) :\n", " intMono[d] = Simpson( a, b, N, monomial )" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# generating random coefficients \n", "coefs = np.random.rand(maxD+1)\n", "\n", "# evaluating Simpson on the polynomial :\n", "intPoly = Simpson(a,b,N, p)\n", "\n", "# computing integral by linearity. Remeber that in polyval, the order of the coefficients is opposite !\n", "intSum = 0\n", "for d in range(maxD+1) :\n", " intSum = intSum + coefs[maxD-d]*intMono[d]\n", "\n", "print(f'Simpson on a_{d} + a_{d-1} x + ... + a_0x^{d} by linearity : \\n'\n", " f'\\t {intPoly:.6f} - {intSum:.6f} = {intPoly-intSum:.6e}')\n", "\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Exercice Série 6, Ex 6 : Degré d'exactitude d'une formule de quadrature" ] }, { "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", "\n", "# In my case, the IntegrationLib is in the parent directory,\n", "# therfore have have to add the aprent directory to path :\n", "import sys \n", "sys.path.append('..')\n", "\n", "from IntegrationLib import *\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 = (b-a) * f( (a+b)/2 )\n", "intTrap = (b-a) * ( f(a)+f(b) )/2\n", "intSimp = (b-a)/6 * ( f(a) + 4*f((a+b)/2) + f(b) )\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", "N = 2**np.linspace(2,7,6).astype(int)\n", "\n", "for n in N :\n", " errmp.append( np.abs( Midpoint( a, b, n, f) - Iexact ) ) \n", " errtrap.append( np.abs( Trapezoidal( a, b, n, f) - Iexact ) ) \n", " errSim.append( np.abs( Simpson( a, b, n, f) - Iexact ) ) \n", " \n", "H = (b-a)/N\n", "plt.plot(H, errmp, 'b:.', H, errtrap, 'c:*', H, errSim, 'g:*')\n", "\n", "plt.plot(H, H**2 * (errmp[0]/H[0]**2)*5, 'k:', H, H**4 * (errSim[0]/H[0]**4)*5, 'k:')\n", "\n", "plt.legend(['rectangle', 'trapeze', 'Simpson', '$H^2$', '$H^4$'])\n", "plt.xlabel('H'); plt.ylabel('err');\n", "\n", "plt.xscale('log') \n", "plt.yscale('log')\n", "plt.grid(True)\n", " " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "*Commentaires :*\n", "\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ - "slopeMP = ( np.log(errmp[-1]) - np.log(errmp[0]) ) / ( np.log(N[-1]) - np.log(H[0]) ) \n", - "slopeTrap = ( np.log(errtrap[-1]) - np.log(errtrap[0]) ) / ( np.log(N[-1]) - np.log(H[0]) ) \n", - "slopeSim = ( np.log(errSim[-1]) - np.log(errSim[0]) ) / ( np.log(N[-1]) - np.log(H[0]) ) \n", + "slopeMP = ( np.log(errmp[-1]) - np.log(errmp[0]) ) / ( np.log(H[-1]) - np.log(H[0]) ) \n", + "slopeTrap = ( np.log(errtrap[-1]) - np.log(errtrap[0]) ) / ( np.log(H[-1]) - np.log(H[0]) ) \n", + "slopeSim = ( np.log(errSim[-1]) - np.log(errSim[0]) ) / ( np.log(H[-1]) - np.log(H[0]) ) \n", "\n", "print(f'La convergence numérique est environ de ')\n", "print(f'Rectangle : {slopeMP:.2f}')\n", "print(f'Trapeze : {slopeTrap:.2f}')\n", "print(f'Simpson : {slopeSim:.2f}')" ] }, { "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", "# recording for which degrees the integral s exact (up to epsilon)\n", "exactDegree = -1\n", "epsilon = 1e-12\n", "\n", "for d in range(7) :\n", " intsim = Midpoint( a, b, N, monomial )\n", " intExact = 1/(d+1)\n", " print(f'Midpoint on x^{d} : {intsim:.6f} - {intExact:.6f} = {intsim-intExact:.6e}')\n", " if np.abs(intsim-intExact) < epsilon :\n", " exactDegree = d\n", " \n", "print(f'Midpoint is exact up to degree {exactDegree}')" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "for d in range(7) :\n", " intsim = Trapezoidal( a, b, N, monomial )\n", " intExact = 1/(d+1)\n", " print(f'Trapezoidal on x^{d} : {intsim:.6f} - {intExact:.6f} = {intsim-intExact:.6e}')\n", " if np.abs(intsim-intExact) < epsilon :\n", " exactDegree = d\n", " \n", "print(f'Trapezoidal is exact up to degree {exactDegree}')" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "for d in range(7) :\n", " intsim = Simpson( a, b, N, monomial )\n", " intExact = 1/(d+1)\n", " print(f'Simpson on x^{d} : {intsim:.6f} - {intExact:.6f} = {intsim-intExact:.6e}')\n", " if np.abs(intsim-intExact) < epsilon :\n", " exactDegree = d\n", " \n", "print(f'Simpson is exact up to degree {exactDegree}')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "*Commentaires :*" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Exercice Série 6, Ex 7 : Convergence pour fonction non-lisse" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "import matplotlib.pyplot as plt\n", "import numpy as np\n", "\n", "from IntegrationLib import Midpoint\n", "from IntegrationLib import Trapezoidal \n", "from IntegrationLib import Simpson\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", "errmp = []\n", "errtrap = []\n", "errSim = []\n", "\n", "N = [1,9,99,999]\n", "#N = [2,10,100,1000]\n", "\n", "for i in range(4) :\n", " errmp.append( np.abs( Midpoint( a, b, N[i], f) - Iexact ) ) \n", " errtrap.append( np.abs( Trapezoidal( a, b, N[i], f) - Iexact ) ) \n", " errSim.append( np.abs( Simpson( a, b, N[i], f) - Iexact ) ) \n", " " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Partie 3" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "H = 2./np.array(N)\n", "plt.plot(H, errmp, 'b:.', H, errtrap, 'c:*', H, errSim, 'g:*')\n", "\n", "plt.plot(H, H**2 * (errmp[0]/H[0]**2)*5, 'k:', H, H**4 * (errSim[0]/H[0]**4)*5, 'k:')\n", "\n", "plt.legend(['rectangle', 'trapeze', 'Simpson', '$H^2$', '$H^4$'])\n", "plt.xlabel('H'); plt.ylabel('err');\n", "\n", "plt.xscale('log') \n", "plt.yscale('log')\n", "plt.grid(True)\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "*Commentaires :*" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Partie 4" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "slopeMP = ( np.log(errmp[-1]) - np.log(errmp[0]) ) / ( np.log(H[-1]) - np.log(H[0]) ) \n", "slopeTrap = ( np.log(errtrap[-1]) - np.log(errtrap[0]) ) / ( np.log(H[-1]) - np.log(H[0]) ) \n", "slopeSim = ( np.log(errSim[-1]) - np.log(errSim[0]) ) / ( np.log(H[-1]) - np.log(H[0]) ) \n", "\n", "print(f'La convergence numérique est environ de ')\n", "print(f'Rectangle : {slopeMP:.2f}')\n", "print(f'Trapeze : {slopeTrap:.2f}')\n", "print(f'Simpson : {slopeSim:.2f}')\n", "\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "*Commentaires :*" ] }, { "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, le point $x=0$ est un point $x_i$ (pour $i=N/2$) exactement entre deux sous-intervalles.\n", "La fonction $f$ est régulière à droite et à gauche, en particulier, on retrouve l'ordre $4$ par rapport à\n", "$H$ de chaque côté.\n", "\n", "* si $N$ est impair, le point $x=0$ est au milieu d'un sous-intervalle.\n", "La fonction $f$ n'est pas régulière dans cet intervalle, ce qui donnt un ordre de convergence plus petit." ] }, { "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 (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.7.16" } }, "nbformat": 4, "nbformat_minor": 4 }