Dérivation et intégration numériques

Objectifs du chapitre

À la fin de ce chapitre, vous serez capable de:

• établir par la formule de Taylor l'ordre d'une formule de différences finies, et distinguer la différence avant, d'ordre 1, de la différence centrée, d'ordre 2;
• calculer le pas $h$ qui minimise la somme de l'erreur de troncature et de l'erreur d'arrondi d'une formule de dérivation, et prédire la précision maximale atteignable — puis la vérifier en exécutant le calcul;
• construire une formule de quadrature comme l'intégrale d'un interpolant, et déterminer son degré d'exactitude;
• démontrer les formules d'erreur du trapèze composite et de Simpson composite, avec leurs constantes, et en déduire le nombre de sous-intervalles nécessaire à une précision donnée;
• accélérer une formule composite par l'extrapolation de Richardson et construire le tableau de Romberg;
• énoncer et démontrer que la quadrature de Gauss à $n$ nœuds est exacte pour les polynômes de degré $2n-1$, appliquer les formules à 2 et 3 nœuds sur un intervalle quelconque, et reconnaître les intégrales — impropres, oscillantes — où ces formules échouent.

Deux opérations, deux destins

Dériver et intégrer sont, en analyse, des opérations symétriques. Numériquement, elles ne le sont pas du tout, et c'est la première chose que ce chapitre doit faire comprendre.

L'intégration est une opération régularisante: elle moyenne. Une perturbation de la fonction d'amplitude $\eta$ change son intégrale sur $[a,b]$ d'au plus $(b-a)\eta$, quelle que soit la forme de la perturbation. Les erreurs de signes opposés se compensent, et une intégrale calculée à partir de valeurs bruitées reste raisonnable.

La dérivation est une opération amplificatrice: elle divise par $h$. Une perturbation d'amplitude $\eta$ sur les valeurs de $f$ change un quotient de différences de $\eta /h$ au moins, et ce facteur $1/h$ explose précisément lorsqu'on cherche à améliorer la formule en raffinant le pas. Le problème «dériver une fonction connue par des valeurs approchées» est mal conditionné, au sens du chapitre 1, et aucun algorithme ne le rendra bien conditionné.

C'est pourquoi les deux moitiés de ce chapitre ne racontent pas la même histoire. Dans la première, l'erreur de troncature et l'erreur d'arrondi se heurtent frontalement, la précision atteignable plafonne bien au-dessus de ${\epsilon }_{\text{mach}}$, et l'essentiel est de savoir où s'arrêter. Dans la seconde, l'erreur d'arrondi reste tapie à $10^{-16}$ pendant que l'erreur de troncature descend de six, huit, douze ordres de grandeur, et l'essentiel est de choisir la formule qui descend le plus vite.

Le fil conducteur numérique du chapitre est l'intégrale témoin du cours,

Cette intégrale n'a pas de primitive élémentaire: la fonction $x\mapsto {\mathrm{e}}^{-x^2}$ est continue, indéfiniment dérivable, bornée, et pourtant aucune combinaison finie de fonctions usuelles n'a ${\mathrm{e}}^{-x^2}$ pour dérivée — c'est un théorème de Liouville, et c'est la raison d'être de ce chapitre. La valeur de référence (8.1) n'a donc pas été obtenue par une primitive, mais par la série entière, intégrée terme à terme sur $[0,1]$:

une série alternée dont les termes décroissent très vite: quarante termes suffisent à saturer la double précision, et tous les chiffres de (8.1) sont ceux que la somme exacte en fractions rationnelles produit. C'est cette valeur qui servira de vérité de référence dans toutes les tables du chapitre.

Dérivation numérique

La différence avant et son ordre 1

On dispose d'une fonction $f$ que l'on sait évaluer, et l'on veut $f^{\prime }(x_0)$ sans calculer la dérivée analytiquement — parce que $f$ est le résultat d'une simulation, d'une table de mesures ou d'un programme dont on n'a pas l'expression. La définition de la dérivée suggère immédiatement un candidat.

Que valent-elles? La réponse tient entièrement dans la formule de Taylor avec reste de Lagrange (Analyse I, chapitre 9), qui est l'unique outil de toute cette section.

Démonstration. La formule de Taylor à l'ordre 1 avec reste de Lagrange s'écrit, pour $f\in {\mathcal{C}}^2$: il existe $\xi \in ]x_0,x_0+h[$ tel que

On soustrait $f(x_0)$ des deux côtés et l'on divise par $h\ne 0$:

ce qui est exactement (8.4). La majoration s'en déduit en remplaçant $|f^{\prime \prime }(\xi )|$ par son maximum sur l'intervalle. $\square$

Le terme $\frac{h}{2}{f}^{\prime \prime }(\xi )$ est l'erreur de troncature de la formule, au sens exact du chapitre 1: c'est ce qu'on perd en remplaçant une limite par l'un de ses termes, et il resterait identique sur une machine arithmétiquement parfaite. Le mot «ordre 1» signifie une chose très concrète: diviser $h$ par dix divise l'erreur par dix, et pas davantage.

La différence centrée et son ordre 2

Prendre le point symétrique à gauche coûte une évaluation de plus et fait gagner un ordre entier. La raison est une annulation, et elle mérite d'être vue.

Démonstration. Écrivons deux développements de Taylor à l'ordre 2 avec reste de Lagrange, l'un vers la droite, l'autre vers la gauche: il existe ${\xi }_{+}\in ]x_0,x_0+h[$ et ${\xi }_{-}\in ]x_0-h,x_0[$ tels que

En soustrayant, les termes en $f(x_0)$ et en $f^{\prime \prime }(x_0)$ disparaissent, parce qu'ils sont portés par des puissances paires de $h$ et entrent donc avec le même signe dans les deux lignes:

La fonction $f^{\prime \prime \prime }$ étant continue, la valeur intermédiaire $\frac{1}{2}(f^{\prime \prime \prime }({\xi }_{+})+f^{\prime \prime \prime }({\xi }_{-}))$, qui est comprise entre et sur , est atteinte en un point de cet intervalle (théorème des valeurs intermédiaires). Donc

et en divisant par $2h$ on obtient (8.5). $\square$

Trois autres formules classiques s'obtiennent de la même manière et seront utilisées plus loin ou dans les exercices:

formule décentrée d'ordre 2, utile au bord d'un intervalle où $f(x_0-h)$ n'existe pas, et

la différence seconde centrée, d'ordre 2, qui est le cœur de toutes les discrétisations d'équations aux dérivées partielles du second ordre.

La collision: troncature contre arrondi

Tout ce qui précède est vrai sur une machine exacte. Sur une vraie machine, les valeurs $f(x_0)$ et $f(x_0+h)$ ne sont pas connues exactement: chacune est entachée d'une erreur relative d'au moins l'unité d'arrondi $u={\epsilon }_{\text{mach}}/2=2^{-53}\approx 1,11\cdot 10^{-16}$, et souvent davantage si est elle-même le résultat d'un calcul. Notons la valeur effectivement obtenue et supposons

au voisinage de $x_0$. L'erreur du quotient de différences calculé se décompose alors en deux morceaux:

Le second est borné par $\frac{M_2}{2}h$, d'après le théorème 8.1. Le premier est borné par $2u{M}_0/h$, puisque le numérateur porte deux erreurs d'au plus $u{M}_0$ chacune et qu'on divise par . D'où la majoration qui gouverne toute la dérivation numérique:

Les deux termes vont en sens opposés. Quand $h$ diminue, le premier décroît linéairement et le second croît comme $1/h$. Leur somme ne peut donc pas tendre vers zéro: elle passe par un minimum, et ce minimum est la précision maximale que la formule peut atteindre, quelle que soit la patience de l'utilisateur.

Démonstration. La fonction $\phi$ est dérivable sur $]0,+\infty [$ et

Elle s'annule en l'unique point positif $h^{\ast }$ tel que $h^{\ast 2}=4u{M}_0/M_2$, c'est-à-dire $h^{\ast }=2\sqrt{u{M}_0/M_2}$. Comme , la fonction est strictement convexe et ce point critique est bien son minimum global. Sa valeur s'obtient en substituant:

Les deux termes sont égaux au point optimal: à $h^{\ast }$, la troncature et l'arrondi contribuent exactement autant. C'est la lecture qualitative à retenir, et elle vaut pour n'importe quel compromis de cette forme. $\square$

Le même calcul appliqué à la différence centrée, dont la troncature est en $h^2{M}_3/6$ et l'arrondi en $u{M}_0/h$ (un seul $h$ au dénominateur de $2h$, et deux erreurs d'arrondi, d'où le facteur 1 au lieu de 2), donne

soit un pas optimal de l'ordre de $10^{-5}$ à $10^{-6}$ et une précision de l'ordre de ${\epsilon }_{\text{mach}}^{2/3}\approx 4\cdot 10^{-11}$: au lieu de la moitié. Une formule d'ordre plus élevé ne se contente donc pas de converger plus vite, elle relève aussi le plancher.

L'expérience centrale du chapitre

Ces formules ne valent rien tant qu'on ne les a pas confrontées à la machine. Faisons-le sur un cas où toutes les constantes sont connues exactement.

Intégration numérique

On passe à l'autre moitié du chapitre, et l'atmosphère change complètement: ici, l'erreur d'arrondi reste à $10^{-16}$ et l'erreur de troncature descend aussi bas qu'on veut, pourvu qu'on choisisse la bonne formule.

Une formule de quadrature est l'intégrale d'un interpolant

Le procédé est toujours le même. On ne sait pas intégrer $f$; on sait intégrer les polynômes. On remplace donc $f$ par un polynôme qui lui ressemble, et l'on intègre celui-là exactement.

L'égalité (8.13) mérite d'être soulignée: les poids ne dépendent pas de $f$. Ils ne dépendent que des nœuds et de l'intervalle. On les calcule une fois, on les tabule, et la formule s'applique ensuite à n'importe quelle fonction. C'est précisément ce que font toutes les tables de quadrature.

Le degré d'exactitude est la mesure de qualité d'une formule, et il se teste par un calcul fini: il suffit de vérifier l'égalité sur les monômes $1,x,x^2,\dots$ jusqu'à ce qu'elle tombe en défaut, puisque l'intégrale et $Q$ sont toutes deux linéaires.

Démonstration. Supposons $Q$ interpolatoire et soit $p$ un polynôme de degré $\le n$. Son polynôme d'interpolation aux $n+1$ nœuds est $p$ lui-même, par unicité de l'interpolant (chapitre 6): deux polynômes de degré $\le n$ qui coïncident en $n+1$ points distincts sont égaux. Donc $Q(p)={\int }_a^{b}p={\int }_a^{b}p$, et le degré d'exactitude est au moins .

Réciproquement, supposons le degré d'exactitude au moins $n$. Les polynômes de Lagrange ${\ell }_i$ sont de degré $n$, donc $Q({\ell }_i)={\int }_a^b{\ell }_i$. Or , d'où . Les poids valent donc bien (8.13) et la formule est interpolatoire.

Rectangles, point milieu, trapèze, Simpson

Les formules élémentaires s'obtiennent en prenant $n=0,1,2$ et les nœuds les plus naturels. Posons $h=b-a$ pour les formules à un et deux nœuds, et $m=(a+b)/2$.

Ces degrés ont été vérifiés numériquement, monôme par monôme, sur $[0,1]$: la formule du point milieu intègre $1$ et $x$ exactement et se trompe sur $x^2$ ($0,25$ au lieu de $1/3$); celle de Simpson intègre exactement $1$, $x$, $x^2$ , et se trompe sur ( au lieu de ).

Deux surprises sont à commenter.

Le point milieu bat le trapèze, alors qu'il coûte une évaluation au lieu de deux. Sur l'intégrale témoin avec un seul intervalle, le trapèze donne $0,6839397206$ (erreur $6,29\cdot 10^{-2}$) et le point milieu $0,7788007831$ (erreur $3,20\cdot 10^{-2}$), soit exactement deux fois mieux. Nous verrons pourquoi: les constantes d'erreur sont et , de signes opposés et dans un rapport 2. Géométriquement, l'excès du trapèze au-dessus d'une courbe concave est le double du défaut du rectangle médian en dessous.

Simpson gagne un degré gratuitement. Son interpolant est une parabole, donc un polynôme de degré 2; le théorème 8.4 ne garantit qu'un degré d'exactitude de 2. Il vaut 3, parce que les nœuds sont symétriques autour de $m$: l'erreur d'interpolation sur $x^3$ est une fonction impaire autour de $m$, et son intégrale sur un intervalle symétrique est nulle. C'est le même mécanisme que celui de la différence centrée, et c'est ce qui fait de Simpson la formule la plus employée de toute l'analyse numérique élémentaire.

Les constantes d'erreur, démontrées

Passons aux formules d'erreur. Elles ne seront pas citées: elles seront obtenues, et la méthode est instructive dans les deux cas.

Démonstration. Soit $p_1$ le polynôme d'interpolation de $f$ en $a$ et $b$. La formule d'erreur d'interpolation (chapitre 6) donne, pour chaque $x\in [a,b]$, l'existence d'un ${\xi }_x\in ]a,b[$ tel que

Comme la formule du trapèze est exactement ${\int }_a^b{p}_1$, on obtient

Le facteur $(x-a)(x-b)$ est de signe constant sur $[a,b]$: il y est négatif ou nul. C'est le point qui permet d'appliquer le théorème de la moyenne généralisé pour les intégrales, à condition que $x\mapsto f^{\prime \prime }({\xi }_x)$ soit continue. Elle l'est: sur , la relation ci-dessus donne

quotient de deux fonctions continues dont le dénominateur ne s'annule pas, et l'on vérifie que ce quotient se prolonge par continuité en $a$ et en $b$ par la règle de l'Hospital. Il existe donc $\xi \in ]a,b[$ tel que

Il ne reste qu'à calculer cette intégrale. Le changement de variable $x=a+t$, $t\in [0,h]$, donne

D'où le résultat: l'erreur vaut $\frac{f^{\prime \prime }(\xi )}{2}\cdot (-\frac{h^3}{6})=-\frac{h^3}{12}{f}^{\prime \prime }(\xi )$.

La formule de Simpson ne se traite pas ainsi, et il faut comprendre pourquoi: son degré d'exactitude est 3, alors que son interpolant n'est que de degré 2; l'erreur d'interpolation ne peut donc pas produire la bonne constante. Il faut un argument qui utilise le degré 3, et la manière la plus élémentaire est de faire varier la longueur de l'intervalle.

Démonstration. Introduisons la fonction qui compare l'intégrale et la formule sur l'intervalle $[m-t,m+t]$, pour $t\in [0,h]$:

On a $F(0)=0$, et $F(h)$ est l'erreur cherchée. Dérivons trois fois, en utilisant le théorème fondamental pour le premier terme.

soit, en regroupant,

Puis

c'est-à-dire

et enfin, les deux premiers termes de la dérivée s'annulant exactement,

La formule des accroissements finis appliquée à $f^{\prime \prime \prime }$ sur $[m-t,m+t]$ (Analyse I, chapitre 8) fournit un point ${\eta }_t$ tel que $f^{\prime \prime \prime }(m+t)-f^{\prime \prime \prime }(m-t)=2t{f}^{(4)}({\eta }_t)$, donc

Notons $\mu ={\min }_{[a,b]}{f}^{(4)}$ et $M={\max }_{[a,b]}{f}^{(4)}$. De (8.16) on tire l'encadrement pour tout . Comme — ce qui se lit directement sur les trois expressions ci-dessus —, on intègre trois fois de à en conservant les inégalités, et comme

il vient

Le nombre $-\frac{90}{h^5}F(h)$ est donc compris entre $\mu$ et $M$; la fonction $f^{(4)}$ étant continue sur $[a,b]$, le théorème des valeurs intermédiaires garantit qu'elle prend cette valeur en un point , ce qui donne (8.15).

Formules composites: de $O(h^{p+1})$ à $O(h^p)$

En pratique, on n'applique jamais ces formules sur $[a,b]$ tout entier. On découpe l'intervalle en $n$ sous-intervalles de longueur $h=(b-a)/n$, avec $x_i=a+ih$, et l'on applique la formule élémentaire sur chacun. C'est la formule composite, et c'est la seule forme utilisable.

Le coût est de $n+1$ évaluations de $f$ dans les deux cas: les valeurs aux nœuds intérieurs sont partagées entre sous-intervalles voisins, et c'est la raison pour laquelle le trapèze composite ne coûte pas $2n$ évaluations.

def trapezes(f, a, b, n):
    """Trapeze composite sur [a, b] avec n sous-intervalles."""
    h = (b - a) / n
    s = 0.5 * (f(a) + f(b))
    for i in range(1, n):
        s += f(a + i * h)
    return h * s


def simpson(f, a, b, n):
    """Simpson composite; n doit etre pair."""
    h 

Les deux fonctions font $n+1$ évaluations de $f$, $n$ additions et une multiplication: le coût est linéaire en $n$, et il est entièrement dominé par les évaluations de $f$ dès que celle-ci est un tant soit peu chère.

Démonstration. Sur chaque sous-intervalle $[x_{i-1},x_i]$, de longueur $h$, le théorème 8.5 fournit un point ${\xi }_i\in ]x_{i-1},x_i[$ tel que

La relation de Chasles donne ${\int }_a^{b}f={\sum }_{i=1}^n{\int }_{x_{i-1}}^{x_i}f$, et la somme des formules élémentaires est exactement , chaque nœud intérieur étant compté deux fois avec le poids . En sommant,

La moyenne $\bar{\mu }$ est comprise entre ${\min }_{[a,b]}{f}^{\prime \prime }$ et ${\max }_{[a,b]}{f}^{\prime \prime }$; comme est continue sur , le théorème des valeurs intermédiaires garantit l'existence d'un tel que . Enfin , donc , ce qui donne (8.19).

Démonstration. Les sous-intervalles se groupent par paires: pour $j=1,\dots ,n/2$, la paire $[x_{2j-2},x_{2j}]$ a pour longueur $2h$, donc pour demi-longueur $h$, et son point milieu est . Le théorème 8.6 s'y applique et fournit tel que

La somme sur $j$ des formules élémentaires reconstitue exactement (8.18): les nœuds d'indice impair reçoivent le coefficient 4 une seule fois, les nœuds d'indice pair intérieurs reçoivent 1 de chaque côté, soit 2, et les extrémités reçoivent 1. D'où

la dernière égalité par le même argument de valeur intermédiaire que ci-dessus, $f^{(4)}$ étant continue. Comme $nh=b-a$,

ce qui est (8.20). $\square$

L'intégrale témoin, ses erreurs et ses rapports

Passons à la mesure. Voici ce que les deux formules donnent sur (8.1), chaque nombre étant celui que le programme a imprimé, et chaque rapport celui qu'il a calculé — pas la valeur théorique.

La colonne des rapports est la preuve. Pour le trapèze, elle vaut $4,02$, $4,01$, $4,00$, $4,00$: doubler $n$ divise $h$ par deux et l'erreur par $2^2=4$, conformément à (8.19). Pour Simpson, elle tend vers 16, c'est-à-dire $2^4$, conformément à (8.20).

Mais elle ne vaut pas 16 tout de suite, et c'est l'information honnête de cette table. Au premier raffinement, le rapport de Simpson vaut 11,40, pas 16. Il faut dire pourquoi. La formule (8.20) s'écrit $-\frac{(b-a)h^4}{180}{f}^{(4)}({\xi }_n)$, où le point ${\xi }_n$ . Le rapport de deux erreurs consécutives vaut donc

et ce second facteur ne vaut 1 qu'à la limite. Pour $f(x)={\mathrm{e}}^{-x^2}$, la dérivée quatrième $f^{(4)}(x)=(16x^4-48x^2+12){\mathrm{e}}^{-x^2}$ varie entre et sur , en changeant deux fois de signe: à grossier, les sautent d'une région à l'autre et le rapport en pâtit. En remontant la formule depuis les erreurs mesurées, on trouve pour , puis pour et pour : la valeur se stabilise, et c'est en même temps que le rapport s'approche de 16. , et arrondir 11,40 en 16 serait dissimuler exactement ce que la table a à enseigner.

Vérifions enfin que les bornes tiennent. Sur $[0,1]$, $\max |f^{\prime \prime }|=2$ (atteint en $x=0$) et $\max |f^{(4)}|=12$ (également en $x=0$). Pour , :

Les deux bornes sont respectées, avec un facteur 3 à 8 de marge — c'est ce qu'on attend d'une majoration qui remplace $f^{(k)}(\xi )$ par son maximum.

Le point milieu composite, et une identité utile

La formule du point milieu composite,

n'utilise aucune valeur aux extrémités, ce qui la rend précieuse pour les intégrandes singulières au bord. Son erreur est

d'ordre 2 comme le trapèze, avec une constante deux fois plus petite et de signe opposé. Sur l'intégrale témoin, le rapport mesuré des deux erreurs vaut $1,988$ pour $n=2$, puis $1,997$, $1,999$ et $2,000$: il tend bien vers 2, et le trapèze sous-estime pendant que le point milieu surestime.

De là découle une identité que l'on retrouvera à la section suivante:

Vérification numérique sur l'intégrale témoin, avec $n=2$: $M_2=0,754597943772$ et $T_2=0,731370251829$, donc , qui est exactement lu dans la table ci-dessus. La combinaison élimine le terme en puisque les deux constantes sont dans le rapport : c'est déjà, sans le nom, de l'extrapolation.

L'extrapolation de Richardson

Les formules d'erreur ont une conséquence pratique remarquable: si l'on connaît la forme de l'erreur, deux calculs approchés en donnent un troisième, bien meilleur, sans aucune évaluation supplémentaire de $f$.

Démonstration. Par hypothèse, $A(h)=A+c{h}^p+O(h^q)$ et $A(h/2)=A+c{h}^p/2^p+O(h^q)$. Multiplions la seconde par et soustrayons la première:

Le terme en $h^p$ s'est annulé exactement, parce que le coefficient $c$ est le même dans les deux développements. Il ne reste qu'à diviser par $2^p-1\ne 0$. $\square$

Appliquée au trapèze composite ($p=2$), l'extrapolation donne

et un petit calcul montre que cette combinaison est exactement $S_{2n}$: l'extrapolation de Richardson du trapèze est la formule de Simpson. C'est la même identité que (8.23), sous un autre habillage, et elle se vérifie sur la table de la page précédente: $(4\times 0,7430-0,7314)/3$ redonne bien la valeur de Simpson.

Pour aller plus loin il faut savoir que le développement de l'erreur du trapèze ne contient que des puissances paires de $h$. C'est le contenu de la formule d'Euler–Maclaurin:

où les $c_k$ dépendent de $f$ et de ses dérivées aux extrémités, mais pas de $h$. Nous l'admettons: sa démonstration passe par les polynômes et les nombres de Bernoulli et relève d'un cours sur les développements asymptotiques. Ce que nous en utilisons est uniquement la structure — des puissances paires — et c'est elle qui rend l'extrapolation si efficace ici: chaque application élimine $h^2$, puis $h^4$, puis $h^6$, et fait donc gagner ordres à chaque fois, pas un.

Le tableau de Romberg

L'application répétée de (8.24) aux trapèzes s'organise en un tableau triangulaire.

Chaque nouvelle ligne ne coûte que $2^{k-1}$ évaluations nouvelles, puisque les nœuds de $T_{2^{k-1}}$ sont la moitié de ceux de $T_{2^k}$. Voici le tableau complet sur l'intégrale témoin, tel que le programme l'a imprimé.

Et voici les erreurs correspondantes, qui sont la vraie lecture du tableau:

Trois remarques.

Le gain est spectaculaire et gratuit. La dernière ligne utilise 17 évaluations de $f$ — exactement celles de $T_{16}$ — et la meilleure valeur du tableau, $R_{4,3}$, vaut $0,746824132647$, à $1,65\cdot 10^{-10}$ de la vérité, soit que le trapèze qui a fourni les mêmes valeurs. Les 16 combinaisons linéaires qui séparent les deux ne coûtent rien.

La première colonne d'extrapolation, c'est Simpson. Comparez $R_{k,1}$ avec la colonne Simpson de la table précédente: $0,747180428910$, $0,746855379791$, $0,746826120527$ — ce sont exactement $S_2$, , , chiffre pour chiffre. Rien de mystérieux: c'est l'identité démontrée plus haut.

Le tableau n'est pas monotone, et c'est normal. À la ligne $n=8$, la colonne 3 ($1,14\cdot 10^{-7}$) est moins bonne que la colonne 2 ($3,71\cdot 10^{-8}$); à $n=16$ l'ordre se rétablit. Il faut s'y attendre: l'extrapolation élimine le terme dominant du développement (8.25) en supposant que les suivants sont négligeables, ce qui n'est vrai qu'asymptotiquement. À grossier, une extrapolation de trop peut amplifier les termes restants. En pratique, on arrête d'extrapoler quand deux éléments consécutifs de la diagonale cessent de se rapprocher — ce qui, sur ce tableau, se produit vers .

La quadrature de Gauss

Toutes les formules précédentes fixent les nœuds d'avance — équidistants — et n'optimisent que les poids. Avec $n$ nœuds et $n$ poids, on dispose pourtant de $2n$ paramètres libres. Si on les utilisait tous, quel degré d'exactitude pourrait-on atteindre?

L'idée: choisir aussi les nœuds

Traitons le cas $n=2$ à la main, sur $[-1,1]$ (on y ramènera tout intervalle par un changement de variable). Cherchons $x_1,x_2,w_1,w_2$ tels que

pour tout polynôme $p$ de degré $\le 3$. Il suffit de l'imposer sur $1,x,x^2,x^3$, ce qui donne quatre équations pour quatre inconnues:

Le système n'est pas linéaire, mais la symétrie le résout: en cherchant $x_2=-x_1$ et $w_1=w_2=w$, les deux équations impaires sont automatiquement satisfaites, la première donne et la troisième , c'est-à-dire . D'où la :

exacte pour tout polynôme de degré $\le 3$ avec deux évaluations, là où Simpson en demande trois. Le même travail sur $1,x^2,x^4$ avec trois nœuds symétriques $-\tau ,0,\tau$ donne $\tau =\sqrt{3/5}$ et les poids , , :

exacte pour tout polynôme de degré $\le 5$. Ces degrés ont été vérifiés monôme par monôme: la formule (8.27) intègre exactement $1,x,x^2,x^3$ et donne $0,194444$ au lieu de $0,2$ sur $x^4$ (ramené à ); la formule (8.28) est exacte jusqu'à et donne au lieu de sur .

Le schéma «$n$ nœuds, degré $2n-1$» n'est pas une coïncidence, et il ne peut pas être dépassé.

Le théorème

Les racines de $L_2$ sont $\pm 1/\sqrt{3}$ et celles de $L_3$ sont $0$ et : ce sont les nœuds trouvés ci-dessus. Le théorème explique pourquoi.

Démonstration. (1) Soient $y_1,\dots ,y_k$ les points de $]-1,1[$ où $L_n$ , et supposons . Posons , polynôme de degré . Le produit ne change alors de signe nulle part sur — chaque changement de signe de est compensé par celui du facteur correspondant de — donc il est de signe constant, et il n'est pas identiquement nul. Par conséquent

ce qui contredit l'orthogonalité de $L_n$ à tout polynôme de degré $<n$. Donc $k=n$: $L_n$ change de signe en $n$ points distincts de $]-1,1[$, et comme il est de degré , ce sont toutes ses racines, et elles sont simples.

(2) Soit $p$ un polynôme de degré $\le 2n-1$. La division euclidienne par $L_n$ s'écrit

D'une part, l'orthogonalité de $L_n$ à $q$ (de degré $<n$) donne

D'autre part, $L_n(x_i)=0$ pour chaque nœud, donc $p(x_i)=q(x_i)L_n(x_i)+r(x_i)=r(x_i)$ et

Enfin, la formule est interpolatoire à $n$ nœuds, donc exacte pour les polynômes de degré $\le n-1$ (théorème 8.4), et $r$ en est un: $G_n(r)={\int }_{-1}^{1}r$. En recollant les trois égalités, .

(3) Soit une formule quelconque à $n$ nœuds $z_1,\dots ,z_n$ et poids $v_1,\dots ,v_n$. Posons

polynôme de degré $2n$. Il s'annule en chaque nœud, donc ${\sum }_i{v}_{i}p(z_i)=0$. Mais $p$ est un carré non identiquement nul et continu, donc ${\int }_{-1}^{1}p>0$. La formule n'est pas exacte sur .

Changement d'intervalle

Les tables de Gauss sont données sur $[-1,1]$ une fois pour toutes. Le transport vers $[a,b]$ se fait par la bijection affine

qui envoie $[-1,1]$ sur $[a,b]$, d'où

Les poids sont donc multipliés par $(b-a)/2$ et les nœuds transportés. Sur $[0,1]$, la formule à 3 nœuds devient: nœuds $0,112701665379$, $0,5$ et $0,887298334621$, poids $5/18$, $4/9$ et , dont la somme vaut exactement 1, c'est-à-dire — un contrôle de recette immédiat et qui attrape la plupart des fautes de transport.

from math import sqrt

NOEUDS_3 = (-sqrt(3.0 / 5.0), 0.0, sqrt(3.0 / 5.0))
POIDS_3 = (5.0 / 9.0, 8.0 / 9.0, 5.0 / 9.0)


def gauss3(f, a, b):
    """Quadrature de Gauss a 3 noeuds sur [a, b]: 3 evaluations de f."""
    c = 0.5 * (a 

Trois évaluations de $f$, trois multiplications et deux additions: le coût est celui de Simpson à un seul intervalle.

Gauss contre Simpson, sur l'intégrale témoin

Voici la comparaison demandée, chaque valeur ayant été calculée.

La lecture demande de la précision, et c'est pour cela que la colonne des évaluations est là.

À nombre d'évaluations égal, Gauss écrase Simpson. Avec les mêmes trois valeurs de $f$, la formule de Gauss est 37 fois plus précise. C'est la traduction directe du théorème 8.10: Simpson à trois nœuds a un degré d'exactitude de 3, Gauss à trois nœuds un degré d'exactitude de 5.

Mais Gauss à 3 nœuds ne bat pas Simpson à 8 sous-intervalles, et il faut le dire: $9,55\cdot 10^{-6}$ contre $1,99\cdot 10^{-6}$, soit un facteur 4,8 en faveur de Simpson. Ce serait étonnant s'il n'utilisait pas trois fois plus d'évaluations. La comparaison honnête est celle qui fixe le budget: à 4 évaluations, Gauss atteint $3,35\cdot 10^{-7}$, c'est-à-dire . À budget égal ou moindre, Gauss gagne toujours; à budget triple, Simpson rattrape. C'est exactement ce qu'on attend de deux méthodes dont l'une a un ordre supérieur.

On peut aussi composer Gauss: appliquer (8.27) sur chacun des $n$ sous-intervalles d'une subdivision. Le résultat est une formule d'ordre 4, comme Simpson, mais avec une bien meilleure constante. Sur l'intégrale témoin, avec $n$ sous-intervalles et $2n$ évaluations:

Le rapport tend vers 16: ordre 4 confirmé, avec la même montée progressive que pour Simpson et pour la même raison. À 8 évaluations, cette formule fait $1,32\cdot 10^{-6}$, là où Simpson en demande 9 pour $1,99\cdot 10^{-6}$.

Quadrature adaptative, en un paragraphe

Toutes les formules composites ci-dessus répartissent les nœuds uniformément, ce qui est absurde quand la fonction est plate sur la moitié de l'intervalle et brutale sur un petit morceau: on paie partout le pas exigé par la région la plus difficile. La quadrature adaptative corrige cela. Sur un intervalle $[a,b]$, on calcule une approximation $Q_1$ par une formule (Simpson, par exemple), puis $Q_2$ en appliquant la même formule séparément sur $[a,m]$ et $[m,b]$. Comme Simpson est d'ordre 4, la différence estime l'erreur de à un facteur près — c'est encore l'argument de Richardson. Si cette estimation est inférieure à la tolérance demandée, on accepte ; sinon on sur les deux moitiés, avec chacune la moitié de la tolérance. L'algorithme concentre ainsi ses évaluations là où la fonction l'exige, et il rend un résultat accompagné d'une estimation de son erreur — deux propriétés qui expliquent que toute bibliothèque sérieuse de quadrature en une dimension (, ) soit adaptative. Sa faiblesse est son critère: l'estimation locale peut être trompeuse si la formule est accidentellement exacte sur le morceau testé, et il faut donc toujours imposer une profondeur de récursion maximale.

Synthèse

• La dérivation numérique amplifie, l'intégration numérique moyenne. Le quotient de différences divise par $h$ une différence de nombres voisins: c'est un problème mal conditionné au sens du chapitre 1. L'intégration, au contraire, ne peut pas amplifier une perturbation de plus du facteur $b-a$.
• Les formules de différences finies s'obtiennent toutes par la formule de Taylor avec reste de Lagrange. La différence avant est d'ordre 1 avec l'erreur $\frac{h}{2}{f}^{\prime \prime }(\xi )$; la différence centrée est d'ordre 2 avec l'erreur $\frac{h^2}{6}{f}^{\prime \prime \prime }(\xi )$, le gain d'un ordre venant de la seule symétrie, qui annule tous les termes pairs du développement.

Exercices

Vous pouvez afficher le corrigé directement sous chaque énoncé après avoir cherché la solution.