{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# 3 Intégration numérique"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## 3.1 Intégration Numérique\n",
""
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## 3.2 Formules de quadrature sur $[-1,1]$\n",
"\n",
"On cherche à approximer l'intégrale $\\displaystyle\\int_{-1}^1 g(t) dt$ d'une fonction continue.\n",
"\n",
" \n",
"Rappel: **Base de Lagrange**\n",
"\n",
"* $n+1$ noeuds $-1\\leq t_0\\leq t_1 \\leq \\cdots \\leq t_n \\leq 1$.\n",
"* $n+1$ fonctions de base associées $\\left\\{\\varphi_0, \\ldots \\varphi_n\\right\\}$\n",
"\n",
"**Interpolation de Lagrange :** $\\displaystyle \\Pi_ng(t)=\\sum_{j=0}^n g(t_j) \\varphi_j(t)$\n",
"\n",
"$$\\int_{-1}^1 g(t) dt \\approx \\int_{-1}^1 \\Pi_ng(t) dt = \\int_{-1}^1 \\sum_{j=0}^n g(t_j) \\varphi_j(t) dt\n",
" = \\sum_{j=0}^n g(t_j) \\underbrace{\\int_{-1}^1 \\varphi_j(t)dt}_{\\omega_j}$$\n",
"\n",
"On pose $M=n+1$ et à l'aide de ces ingrédients on définit les formules de quadrature.\n",
"\n",
"Une **formule de quadrature** $J(\\cdot)$ permet d'approcher $\\int_{-1}^1 g(t) dt$\n",
"pour une fonction continue $g:[-1,1]\\rightarrow\\rr$. Étant donné $M>0$,\n",
"\n",
"* $M$ points d'intégration $-1\\leq t_1 \\leq \\cdots \\leq t_M \\leq 1$,\n",
"* $M$ poids $\\omega_1,\\ldots,\\omega_M$,\n",
"\n",
"elle s'écrit\n",
"\n",
"$$J(g) = \\sum_{j=0}^n \\omega_j g(t_j).$$\n",
"\n",
"Une formule de quadrature est linéaire (exercice) :\n",
"$$J(f+\\lambda g)= J(f)+\\lambda J(g)$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## 3.2 Quadrature\n",
""
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Degré d'exactitude}\n",
"\n",
"Une formule de quadrature $J$ est **exacte de degré $r$** si pour tout polynome $p$\n",
"de degré $\\leq r$ on a $\\int_{-1}^1 p(t) dt = J(p)$.\n",
"\n",
"\n",
"**Théorème (Poids de quadrature)** Soit $J$ une formule de quadrature avec $M$ noeuds.\n",
"\n",
" $J$ est exacte de degré $M-1$ $\\Longleftrightarrow$ $\\omega_j = \\int_{-1}^1 \\varphi_j(t)dt, \\; j=1,\\ldots, M$,\n",
" où $\\{\\varphi_j, j=1,\\ldots,M\\}$ est la base de Lagrange associée aux noeuds de quadrature.\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Formule de quadrature composite sur $[a,b]$\n",
"\n",
"On cherche à approximer l'intégrale $\\int_{a}^b f(x) dt$ d'une fonction continue sur $[a,b]$\n",
"en utilisant\n",
"\n",
"* une partition en $N$ sous-intervalles de taille $H = \\frac{b-a}N$. Avec\n",
"$x_k = a+kH, k=0,\\ldots,N$ on obtient\n",
"$$\\int_a^b f(x) dx = \\int_{x_0}^{x_1} f(x) dx + \\ldots + \\int_{x_{N-1}}^{x_N} f(x) dx\n",
" = \\sum_{k=0}^{N-1} \\int_{x_{k}}^{x_{k+1}} f(x) dx$$\n",
"\n",
"* le changement de variable\n",
"$$\\int_{x_{k}}^{x_{k+1}} f(x) dx =\n",
" \\int_{-1}^{1} f\\left( \\frac{x_{k}+x_{k+1}}2 + t\\frac{x_{k+1}-x_{k}}2\\right) \\frac{x_{k+1}-x_{k}}2 dt =$$\n",
"\n",
"$$\\int_{-1}^{1} f\\left( \\frac{x_{k}+x_{k+1}}2 + \\frac{H}2t\\right) \\frac{H}2 dt =\n",
" \\frac{H}2 \\int_{-1}^{1} f\\left( x_{k} + \\frac{H}2(t+1)\\right) dt$$\n",
"\n",
"La **formule de quadrature composite** associée à la formule de quadrature $J$ est définie par\n",
"$$L_H(f) = \\frac{H}2 \\sum_{k=0}^{N-1} J(g_k) \\quad\\text{où } g_k(t) = f\\left( x_{k} + \\frac{H}2(t+1)\\right).$$\n",
"\n",
"Par conséquent on a que\n",
"$$L_H(f) = \\frac{H}2 \\sum_{k=0}^{N-1} \\sum_{j=1}^M g_k(t_j)\\omega_j =\n",
" \\frac{H}2 \\sum_{k=0}^{N-1} \\sum_{j=1}^M \\omega_j \\,f\\left( x_{k} + \\frac{H}2(t_j+1)\\right).$$\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## 3.3 Quadrature Composite\n",
""
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## 3.4 Ordre de Convergence\n",
""
]
},
{
"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 4, Série 6\n",
"\n",
"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",
" # size of the subintervals\n",
" H = (b - a) / N\n",
" \n",
" # COMPLETEZ LA SUITE\n",
" \n",
" # quadrature nodes\n",
" # x = \n",
"\n",
" # approximate integral \n",
" # integral = \n",
" return integral\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": [
"(Le résultat précédent doit être 0.3325)\n",
"\n",
"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",
" intmp = midpoint( a, b, N[i], f)\n",
" errmp.append( np.abs( intmp - 1.0/3 ) ) \n",
"\n",
"plt.plot(N, errmp, 'b:.')\n",
"\n",
"plt.xlabel('M'); plt.ylabel('err');\n",
"\n",
"plt.xscale('log') \n",
"plt.yscale('log')\n",
"plt.grid(True)\n",
"plt.title('Le graphe loglog est une droite avec pente égale à -2:')\n",
"plt.show()\n",
"print('Figure: Erreur relative à la méthode du rectangle en échelle logarithmique.\\n\\n')\n",
"\n",
"slope = ( np.log(errmp[4]) - np.log(errmp[0]) ) / ( np.log(N[4]) - np.log(N[0]) ) \n",
"## Essayez de comprendre la signe précédente !\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",
" # size of the subintervals\n",
" H = (b - a) / N\n",
"\n",
" # COMPLETEZ LA SUITE\n",
" \n",
" # quadrature nodes\n",
" # x = \n",
"\n",
" # approximate integral \n",
" # integral = \n",
" return integral\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"errtrap = []\n",
"for i in range(5) :\n",
" inttrap = trapezoidal( a, b, N[i], f)\n",
" errtrap.append( np.abs( inttrap - 1.0/3 ) ) \n",
"\n",
"plt.plot(N, errmp, 'b:.')\n",
"plt.plot(N, errtrap, 'c:*')\n",
"\n",
"\n",
"plt.legend(['rectangle', 'trapeze'])\n",
"plt.xlabel('N'); plt.ylabel('err');\n",
"\n",
"plt.xscale('log') \n",
"plt.yscale('log')\n",
"plt.grid(True)\n",
"plt.title('Les graphes loglog sont des droites avec pente égale à -2:')\n",
"plt.show()\n",
"print('Figure: Erreur relative aux méthodes du rectangle et du trapèze en échelle logarithmique.\\n\\n')\n",
"\n",
"slope = ( np.log(errtrap[4]) - np.log(errtrap[0]) ) / ( np.log(N[4]) - np.log(N[0]) ) \n",
"\n",
"print(f'La convergence numérique est d\\'environ {slope:.2f}')\n",
"\n",
"print(f'En moyenne, le rapport entre les deux erreurs est de {np.mean(np.abs(np.array(errtrap)/np.array(errmp))):.6f}')\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": []
}
],
"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
}