Calcul numérique des intégrales

Comment fait-on pour calculer numériquement une intégrale sur ordinateur ? Deux cas se présentent, selon que l'on ne connaît que certaines valeurs de la fonction (typiquement, par une table de données expérimentales), ou bien que l'on peut l'évaluer en n'importe quel point. Comme nous allons le voir, la deuxième situation est beaucoup plus favorable. Supposons que l'on dispose de $n+1$ abscisses $x_0,\ldots,x_n$ et de $n+1$ ordonnées $y_0,\ldots,y_n$. L'idée la plus naïve est de calculer la somme des aires de rectangles basés sur les intervalles $[x_i,x_{i+1}]$, et de hauteur $y_i$ ou bien $y_{i+1}$ : ce sont les méthodes des rectangles, à gauche, ou à droite.

$$R_g = \sum_{i=1}^ny_{i-1}(x_i-x_{i-1})$$   et$$\quad R_d = \sum_{i=1}^ny_{i}(x_i-x_{i-1})\;.$$

Si on ne sait absolument rien de la fonction à intégrer, il n'y a pas de raison d'aller plus loin. Mais si on imagine un modèle dans lequel les ordonnées $y_i$ sont les évaluations en $x_i$ d'une fonction continue $f$, alors on obtient une meilleure précision par la méthode des trapèzes, qui consiste simplement à prendre la demi-somme de $R_g$ et $R_d$.

$$T = \frac{1}{2}\sum_{i=1}^n(y_{i-1}+y_{i})(x_{i}-x_{i-1})\;.$$

Si on imagine un modèle où la fonction à intégrer $f$ est encore plus lisse, on gagnera en l'interpolant par une parabole sur des triplets de points successifs. C'est la méthode de Simpson.

Pour donner une idée de la précision atteinte, supposons que nous disposions d'une table de $101$ valeurs régulièrement espacées de la fonction sinus sur l'intervalle $[0,\pi/2]$. L'intégrale de $\sin(x)$ sur $[0,\pi/2]$ vaut $1$. Voici ce que donnent la méthode des rectangles à droite et à gauche, la méthode des trapèzes, et la méthode de Simpson.

$$R_g=0.9921255\;,\quad R_d=1.0078334\;,\quad T=0.9999794\;,\quad S=1.0000000003\;.$$

Quand une fonction est donnée par une expression qui permet de la calculer en n'importe quel point, il serait particulièrement inefficace de l'évaluer en $101$ points régulièrement espacés, pour se ramener au cas précédent. On utilise plutôt les méthodes de quadrature de Gauss. Voici comment fonctionne la plus simple, celle de Gauss-Legendre. Soit une fonction $f$ à intégrer sur l'intervalle $[a,b]$. Quitte à effectuer un changement de variable affine, on peut se ramener au cas où $a=-1$ et $b=1$. Les points auxquels on doit évaluer la fonction sont les racines des polynômes de Legendre. On peut définir ces polynômes par récurrence, par une équation différentielle, ou bien comme la dérivée $n$-ième d'un polynôme de degré $2n$.

$$P_n(x)=\frac{1}{2^n(n!)}\Big((x^2-1)^n\Big)^{(n)}\;.$$

Les racines de $P_n$ sont dans l'intervalle $[-1,1]$, et réparties de façon symétrique par rapport à l'origine. Comme $P_n$ est de degré $n$, il a $n$ racines : notons-les $x_1,\ldots,x_n$. À la racine $x_i$, on associe le «poids»$\omega_i$ :

$$\omega_i = \frac{2}{(1-x_i)^2(P'_n(x_i))^2}\;.$$

On calcule ensuite la somme :

$$L_n = \sum_{i=1}^n \omega_i f(x_i)\;.$$

Les racines des polynômes de Legendre, ainsi que les poids qui leur sont associés sont connus depuis longtemps, et inclus dans les bibliothèques de codes des langages de calcul scientifique : le calcul de $L_n$ est donc extrêmement rapide. Vous apprendrez plus tard les raisons mathématiques pour lesquelles ce calcul donne une valeur approchée précise de l'intégrale de $f$. Les résultats sont spectaculaires. Voici pour la même fonction sinus entre 0 et $\pi/2$, ce que donne la méthode de Gauss-Legendre pour $n\leqslant 5$.

$$L_2-1=-1.5 10^{-3}\;,\quad L_3-1=8.1 10^{-6}\;,\quad L_4-1=-2.3 10^{-8}\;,\quad L_5-1=3.9 10^{-11}\;.$$

Ainsi, en évaluant la fonction sinus en $5$ points seulement, on calcule son intégrale avec une précision de l'ordre de $10^{-11}$ : impressionnant non ?