Méthodes de quadrature de Gauss

Pour n donné :

• Les n nœuds x_i sont réels, distincts, uniques et sont les racines du polynôme unitaire de degré n, orthogonal au sous-espace des polynômes de degré n-1 pour le produit scalaire ${\displaystyle \left\langle h,g\right\rangle =\int _{a}^{b}h(x)g(x)\varpi (x)\,\mathrm {d} x}$[2].
• Pour tout i, ${\displaystyle \omega _{i}}$ est égal à ${\displaystyle \int _{a}^{b}l_{i}(x)\varpi (x)\,\mathrm {d} x}$, où l_i est le polynôme interpolateur de Lagrange de degré n-1, prenant la valeur 1 en x_i, et 0 en les x_k pour k différent de i.
• Pour les polynômes de degré inférieur ou égal à 2n-1, il y a égalité entre ${\displaystyle \int _{a}^{b}f(x)\varpi (x)\,\mathrm {d} x}$ et ${\displaystyle \sum _{i=1}^{n}\omega _{i}f(x_{i})}$.

Le domaine d'intégration et la fonction de pondération déterminent le type de la quadrature de Gauss. Le tableau suivant résume les situations les plus communes.

Principales configurations de quadrature de Gauss

Domaine d'intégration (a,b)	Fonction de pondération ${\displaystyle \varpi (x)}$	Famille de polynômes orthogonaux
[–1, 1]	1	Legendre
]–1, 1[	${\displaystyle (1-x)^{\alpha }(1+x)^{\beta }\ ,\ \alpha ,\beta >-1}$	Jacobi
]–1, 1[	${\displaystyle {\frac {1}{\sqrt {1-x^{2}}}}}$	Tchebychev (premier type)
]–1, 1[	${\displaystyle {\sqrt {1-x^{2}}}}$	Tchebychev (second type)
ℝ+	${\displaystyle \mathrm {e} ^{-x}}$	Laguerre
ℝ	${\displaystyle \mathrm {e} ^{-x^{2}}}$	Hermite

Une fois le type de quadrature choisi, la formule à n points s'écrit :

${\displaystyle I(f)=\sum _{i=1}^{n}\omega _{i}f(x_{i}).}$

On définit l'erreur comme ${\displaystyle E(f)=|I-I(f)|}$. Le degré d'exactitude d'une formule de quadrature est le degré le plus élevé de la famille des polynômes annulant E(f). On a le résultat suivant : une formule à n points admet un degré d'exactitude de 2n–1.

Le domaine d'intégration [a, b] peut être changé (au moyen d'un changement de variable) en [–1, 1] avant d'appliquer les méthodes de quadrature de Gauss. Le changement se déroule ainsi :

${\displaystyle \int _{a}^{b}f(x)\,\mathrm {d} x={\frac {b-a}{2}}\int _{-1}^{1}f\left({\frac {b-a}{2}}x+{\frac {a+b}{2}}\right)\,\mathrm {d} x.}$

L'approximation de la valeur de l'intégrale devient :

${\displaystyle {\frac {b-a}{2}}\sum _{i=1}^{n}\omega _{i}f\left({\frac {b-a}{2}}x_{i}+{\frac {a+b}{2}}\right)~.}$

où les x_i sont ici relatifs à l'intervalle [–1, 1].

Pour le problème d'intégration le plus classique, on utilise la méthode de Gauss-Legendre[3]. Il s'agit d'intégrer la fonction f sur le segment [–1, 1]. Les n nœuds sont les racines du n-ième polynôme de Legendre, P_n(x), et les coefficients sont donnés par l'une ou l'autre égalité :

${\displaystyle \omega _{i}={\frac {-2}{(n+1)P'_{n}(x_{i})P_{n+1}(x_{i})}}={\frac {2}{(1-x_{i}^{2})P'_{n}(x_{i})^{2}}}.}$

On peut aussi remarquer que la somme des coefficients est égale à 2. Le tableau suivant donne l'ensemble des informations pour réaliser le calcul approché de I pour les formules à un, deux et trois points.

Nombre de points, n	Poids ${\displaystyle (\omega _{i})}$	Points ${\displaystyle (x_{i})}$	Polynôme de Legendre
1	2	0	x
2	1, 1	–√1/3, √1/3	(3x2 – 1)/2
3	5/9, 8/9, 5/9	–√3/5, 0, √3/5	(5x3 – 3x)/2

On cherche à déterminer ${\displaystyle \int _{-1}^{1}(x+1)^{2}\mathrm {d} x}$. On cherche à intégrer un polynôme de degré 2, 2 points suffisent pour obtenir la valeur exacte.

${\displaystyle \int _{-1}^{1}(x+1)^{2}\mathrm {d} x=1\left({\frac {1}{\sqrt {3}}}+1\right)^{2}+1\left(-{\frac {1}{\sqrt {3}}}+1\right)^{2}={\frac {8}{3}}.}$

On peut facilement vérifier ce résultat car dans cet exemple, on connaît une primitive de (x + 1)² :

${\displaystyle \int _{-1}^{1}(x+1)^{2}\mathrm {d} x=\left[{\frac {(x+1)^{3}}{3}}\right]_{-1}^{1}={\frac {8}{3}}.}$

Cet exemple ne représente pas un cas pratique. En règle générale, on n'obtient jamais un résultat exact et bien entendu, on n'applique pas ces méthodes pour les fonctions dont on connaît une primitive.

Cette formule est associée au poids ${\displaystyle \varpi (x)=(1-x^{2})^{-1/2}}$ sur ]–1, 1[. Pour une formule à n points[4], les nœuds sont

${\displaystyle x_{i}=\cos \left({\frac {(2i-1)\pi }{2n}}\right)}$

et les coefficients :

${\displaystyle w_{i}={\frac {\pi }{n}}.}$

Cette formule est associée au poids ${\displaystyle \varpi (x)={\rm {e}}^{-x}}$ sur ]0, +∞[. Les n nœuds sont les n racines du n-ième polynôme de Laguerre L_n, et les coefficients sont

${\displaystyle w_{i}={\frac {1}{(n+1)L'_{n}(x_{i})L_{n+1}(x_{i})}}.}$

Les coefficients et les nœuds ne peuvent être calculés analytiquement que pour n petit[5]. Par exemple, pour n = 2 :

n	${\displaystyle x_{i}\,}$	${\displaystyle w_{i}\,}$
2	${\displaystyle 2\pm {\sqrt {2}}}$	${\displaystyle {\frac {2\mp {\sqrt {2}}}{4}}}$

Maintenant, pour intégrer une fonction f sur ℝ⁺, il faut remarquer que

${\displaystyle \int _{0}^{+\infty }f(x)\,\mathrm {d} x=\int _{0}^{+\infty }{\frac {f(x)}{\varpi (x)}}\varpi (x)\,\mathrm {d} x.}$

Il reste alors à appliquer la formule de quadrature à la fonction ${\displaystyle g(x)=f(x)/\varpi (x).}$

Sur ℝ, la formule de Gauss-Hermite est caractérisée par la pondération ${\displaystyle \varpi (x)={\rm {e}}^{-x^{2}}}$. Pour une formule à n points[6], les x_i sont calculés comme les n racines du n-ième polynôme d'Hermite H_n ; quant aux pondérations, elles sont obtenues à partir de

${\displaystyle w_{i}={\frac {2^{n+1}n!{\sqrt {\pi }}}{[H_{n}'(x_{i})]^{2}}}.}$

Concernant l'intégration de f sur ℝ, il suffit d'appliquer la formule de quadrature à la fonction ${\displaystyle f(x){\rm {e}}^{x^{2}}.}$

Pour le cas des méthodes de quadrature de Gauss-Lobatto (en) sur l'intervalle [a, b], on impose parmi les r + 1 points de quadrature les deux extrémités de l'intervalle :

${\displaystyle a=x_{1}<x_{2}<\cdots <x_{r+1}=b.}$

Pour un ordre de quadrature r, les points intérieurs de quadrature deviennent alors les zéros du polynôme dérivé du r-1^e polynôme orthogonal :

${\displaystyle \forall i=2,\ldots ,r\ ,\ P_{r-1}'(x_{i})=0.}$

Pour le cas des méthodes de quadrature de Gauss-Radau sur l'intervalle [a, b], on impose parmi les r+1 points de quadrature une des extrémités :

${\displaystyle a=x_{1}<x_{2}<\cdots <x_{r+1}.}$

Par symétrie, on peut également fixer b comme point.

Pour un ordre de quadrature r, les points intérieurs de quadrature deviennent alors les zéros du polynôme :

${\displaystyle {\frac {P_{r-1}(x)+P_{r}(x)}{x-a}}.}$