Problème : la quadrature de Gauss–Legendre

Ce qu'on a le droit d'utiliser. Gram-Schmidt, la meilleure approximation, la base de Lagrange (chapitre 1), la division euclidienne, la factorisation d'un polynôme réel. On ne se servira ni de la formule de Rodrigues ni du théorème de Rolle : l'orthogonalité seule localise les racines.

1) Le polynôme orthogonal unitaire. Existence. Le sous-espace ${\mathbb{R}}_{n-1}[X]$ est de dimension finie ; notons $p$ la projection orthogonale sur lui, et posons $L_n=X^n-p(X^n)$. C'est un polynôme unitaire de degré $n$ (on retranche à $X^n$ un polynôme de degré au plus $n-1$), orthogonal à ${\mathbb{R}}_{n-1}[X]$ par définition de la projection : c'est, à normalisation près, le vecteur d'indice $n$ que produit Gram-Schmidt appliqué à $(1,X,\dots ,X^n)$. Unicité. Si $L$ et $L^{\prime }$ conviennent, $L-L^{\prime }$ est de degré au plus $n-1$ (les termes en $X^n$ se compensent) et orthogonal à ${\mathbb{R}}_{n-1}[X]$ : il est orthogonal à lui-même, donc nul. La norme. Par le théorème de la meilleure approximation, $\Vert L_n\Vert =\Vert X^n-p(X^n)\Vert =d(X^n,{\mathbb{R}}_{n-1}[X])$ : $L_n$ est le polynôme unitaire de degré $n$ de plus petite norme, puisque tout polynôme unitaire de degré $n$ s'écrit $X^n-R$ avec $R\in {\mathbb{R}}_{n-1}[X]$.

Les trois premiers. $L_n$ est le polynôme $P_n$ de l'exercice 5 pour le poids $w=1$ — même caractérisation, donc même polynôme : $L_1=X$, $L_2=X^2-\frac{1}{3}$, $L_3=X^3-\frac{3}{5}X$, de normes au carré $\frac{2}{3}$, $\frac{8}{45}$ et $\frac{8}{175}$, avec $\Vert 1{\Vert }^2=2$. Directement : $⟨X^3,1⟩=⟨X^3,L_2⟩=0$ par imparité, et $\frac{⟨X^3,X⟩}{\Vert X{\Vert }^2}=\frac{2/5}{2/3}=\frac{3}{5}$, donc la projection de $X^3$ sur ${\mathbb{R}}_2[X]$, de base orthogonale $(1,X,L_2)$, est $\frac{3}{5}X$. En normalisant $(1,X,L_2)$, on obtient $e_0=\frac{1}{\sqrt{2}}$, $e_1=\sqrt{\frac{3}{2}}X$ et $e_2=\sqrt{\frac{45}{8}}(X^2-\frac{1}{3})=\frac{\sqrt{10}}{4}(3X^2-1)$.

La distance de $X^3$. $d(X^3,{\mathbb{R}}_2[X]{)}^2=\Vert L_3{\Vert }^2=\frac{8}{175}$, soit $d=\frac{2\sqrt{14}}{35}\approx 0,2138$. Contrôle par Pythagore : $\Vert X^3{\Vert }^2=\frac{2}{7}$, $\Vert \frac{3}{5}X{\Vert }^2=\frac{9}{25}\times \frac{2}{3}=\frac{6}{25}$, et $\frac{2}{7}-\frac{6}{25}=\frac{50-42}{175}=\frac{8}{175}$.

2) Les racines sont réelles, simples, dans l'intervalle. Factorisons $L_n$ dans $\mathbb{R}[X]$ en isolant ses racines situées dans $]-1,1[$ : $$L_n=\prod _r(X-r{)}^{m_r}\times R,$$ où $r$ parcourt les racines de $L_n$ dans $]-1,1[$, de multiplicités $m_r$, et où $R$ n'a aucune racine dans $]-1,1[$. La fonction polynomiale $R$, continue et sans zéro sur $]-1,1[$, y garde un signe constant par le théorème des valeurs intermédiaires. Dans le produit $L_{n}Q$, chaque racine $r$ apparaît avec l'exposant $m_r+1$ si $m_r$ est impair ($r$ est alors l'un des $x_i$) et $m_r$ sinon : tous les exposants sont pairs. Le produit des facteurs $(t-r{)}^{\text{pair}}$ est positif, et $L_{n}Q$ a, sur $]-1,1[$, le signe constant de $R$, puis sur $[-1,1]$ par continuité.

Le polynôme $L_{n}Q$ n'est pas nul : il ne s'annule qu'en un nombre fini de points, et ${\int }_{-1}^1{L}_{n}Q\ne 0$, car une fonction continue de signe constant et d'intégrale nulle serait nulle. Mais si $k\le n-1$, alors $Q\in {\mathbb{R}}_{n-1}[X]$ et ${\int }_{-1}^1{L}_{n}Q=⟨L_n,Q⟩=0$ par le 1) : contradiction. Donc $k\ge n$. Les $x_i$ étant des racines distinctes de $L_n$, qui est de degré $n$, $k\le n$. Donc $k=n$ : $L_n$ a $n$ racines distinctes dans $]-1,1[$, et la somme de leurs multiplicités valant au plus $n$, ses $n$ racines sont réelles, simples, et toutes dans $]-1,1[$.

Le point délicat. L'argument ne demande rien sur la forme de $L_n$ : seulement qu'il soit orthogonal à tout polynôme de degré inférieur. Le polynôme $Q$ « épouse » les changements de signe de $L_n$ ; s'il y en avait moins de $n$, $Q$ serait de degré trop petit pour que $L_n$ puisse lui être orthogonal.

3) La formule à $n$ points. Les $x_i$ sont $n$ réels distincts : les polynômes de Lagrange ${\ell }_i={\prod }_{j\ne i}\frac{X-x_j}{x_i-x_j}$ forment une base de ${\mathbb{R}}_{n-1}[X]$, et tout $P\in {\mathbb{R}}_{n-1}[X]$ s'écrit $P={\sum }_{i=1}^{n}P(x_i){\ell }_i$. En intégrant, par linéarité, $${\int }_{-1}^{1}P=\sum _{i=1}^{n}P(x_i){\int }_{-1}^1{\ell }_i=\sum _{i=1}^n{w}_{i}P(x_i).$$ Rien de remarquable encore : n'importe quels $n$ points distincts donnent, par ce procédé, une formule exacte sur ${\mathbb{R}}_{n-1}[X]$.

4) Le degré double. Soit $P\in {\mathbb{R}}_{2n-1}[X]$. La division euclidienne par $L_n$, de degré $n$, donne $P=L_{n}S+R$ avec $\mathrm{deg}R\le n-1$ et $S\in {\mathbb{R}}_{n-1}[X]$ (car $\mathrm{deg}S=\mathrm{deg}P-n\le n-1$ si $\mathrm{deg}P\ge n$, et $S=0$ sinon). D'une part, par le 1) puis le 3), $${\int }_{-1}^{1}P=⟨L_n,S⟩+{\int }_{-1}^{1}R=0+\sum _{i=1}^n{w}_{i}R(x_i).$$ D'autre part, les $x_i$ étant racines de $L_n$, $P(x_i)=L_n(x_i)S(x_i)+R(x_i)=R(x_i)$. Donc ${\int }_{-1}^{1}P={\sum }_i{w}_{i}P(x_i)$ : la formule est exacte sur ${\mathbb{R}}_{2n-1}[X]$.

Pas au-delà. Soient $y_1,\dots ,y_n$ des réels et ${\alpha }_1,\dots ,{\alpha }_n$ des poids quelconques. Le polynôme $T={\prod }_{i=1}^n(X-y_i{)}^2$, de degré $2n$, est positif et non nul, donc ${\int }_{-1}^{1}T>0$ ; mais la formule lui attribue ${\sum }_i{\alpha }_{i}T(y_i)=0$. Aucune formule à $n$ points n'est exacte sur ${\mathbb{R}}_{2n}[X]$ : le degré $2n-1$ est le meilleur possible, et Gauss l'atteint.

Le point délicat. Tout tient en une ligne : le quotient $S$ est de degré trop petit pour « voir » $L_n$, par orthogonalité, et le reste $R$ est intégré exactement parce qu'il est de petit degré. Placer les nœuds aux racines de $L_n$, c'est ce qui permet d'ignorer le quotient.

5) Des poids positifs. Le polynôme ${\ell }_i^2$ est de degré $2n-2\le 2n-1$ : par le 4), ${\int }_{-1}^1{\ell }_i^2={\sum }_j{w}_j{\ell }_i(x_j{)}^2=w_i$, puisque ${\ell }_i(x_j)$ vaut $1$ si $j=i$ et $0$ sinon. Comme ${\ell }_i^2$ est positif et non nul, $w_i>0$. Enfin, la formule appliquée au polynôme constant $1$ donne ${\sum }_i{w}_i={\int }_{-1}^{1}1=2$. Pourquoi c'est important : si les valeurs de $f$ sont entachées d'erreurs de valeur absolue au plus $\eta$, l'erreur sur le résultat est au plus ${\sum }_i{w}_i\eta =2\eta$ — des poids positifs de somme $2$ n'amplifient pas les erreurs.

6) Deux et trois points. Pour $n=2$, $L_2=X^2-\frac{1}{3}$ a pour racines $x_1=-\frac{1}{\sqrt{3}}$ et $x_2=\frac{1}{\sqrt{3}}$ ; ${\ell }_1=-\frac{\sqrt{3}}{2}(X-\frac{1}{\sqrt{3}})$, d'où $w_1=-\frac{\sqrt{3}}{2}\times (-\frac{2}{\sqrt{3}})=1$, puis $w_2=2-w_1=1$. La formule $P(-\frac{1}{\sqrt{3}})+P\left(\frac{1}{\sqrt{3}}\right)$ est exacte sur ${\mathbb{R}}_3[X]$ — pour $X^2$ : $\frac{1}{3}+\frac{1}{3}=\frac{2}{3}={\int }_{-1}^1{t}^2$ — et fausse pour $X^4$ : $\frac{2}{9}$ au lieu de $\frac{2}{5}$.

Pour $n=3$, $L_3=X(X^2-\frac{3}{5})$ a pour racines $-\sqrt{3/5}$, $0$ et $\sqrt{3/5}$. L'exactitude sur $1$, $X$ et $X^2$, garantie par le 3), donne $w_1+w_2+w_3=2$, $\sqrt{3/5}(w_3-w_1)=0$ et $\frac{3}{5}(w_1+w_3)=\frac{2}{3}$ : d'où $w_1=w_3=\frac{5}{9}$ et $w_2=2-\frac{10}{9}=\frac{8}{9}$. La formule est exacte sur ${\mathbb{R}}_5[X]$ — pour $X^4$ : $2\times \frac{5}{9}\times \frac{9}{25}=\frac{2}{5}$ — et donne $\frac{6}{25}$ au lieu de $\frac{2}{7}$ pour $X^6$.

L'exponentielle. La valeur exacte est $e-e^{-1}\approx 2,350402$.

• Gauss à deux points : $e^{-1/\sqrt{3}}+e^{1/\sqrt{3}}\approx 2,342696$, erreur d'environ $7,7\times 10^{-3}$ ;
• Gauss à trois points : $\frac{5}{9}(e^{-\sqrt{3/5}}+e^{\sqrt{3/5}})+\frac{8}{9}\approx 2,350337$, erreur d'environ $6,5\times 10^{-5}$ ;
• Simpson : $\frac{1}{3}(e^{-1}+4+e)\approx 2,362054$, erreur d'environ $1,2\times 10^{-2}$.

À trois évaluations chacune, Gauss est environ cent quatre-vingts fois plus précis que Simpson ; avec deux seulement, il fait déjà mieux que Simpson avec trois — il atteint le degré $3$ avec deux points, le degré $5$ avec trois.

Le théorème obtenu — la quadrature de Gauss–Legendre (Gauss, 1814). Soit $L_n$ le polynôme unitaire de degré $n$ orthogonal à ${\mathbb{R}}_{n-1}[X]$ pour le produit scalaire ${\int }_{-1}^{1}PQ$. Ses $n$ racines $x_1<\cdots <x_n$ sont réelles, simples et situées dans $]-1,1[$ ; avec les poids $w_i={\int }_{-1}^1{\ell }_i$, strictement positifs et de somme $2$, la formule ${\int }_{-1}^{1}P={\sum }_{i=1}^n{w}_{i}P(x_i)$ est exacte pour tout polynôme de degré au plus $2n-1$, et aucune formule à $n$ points n'est exacte pour tous les polynômes de degré $2n$.

Ce que le problème installe. L'orthogonalité fabrique une méthode numérique : placer les nœuds aux racines d'un polynôme orthogonal double le degré d'exactitude, le maximum possible, avec des poids positifs qui rendent la méthode stable. Deux gestes resserviront : localiser des racines par l'orthogonalité et diviser par le polynôme orthogonal. C'est la méthode d'intégration des logiciels de calcul scientifique, des éléments finis en mécanique et de la chimie quantique.