Extrapolation de Richardson

En analyse numérique, le procédé d'extrapolation de Richardson est une technique d'accélération de la convergence. Il est ainsi dénommé en l'honneur de Lewis Fry RichardsonModèle:,, qui l'a popularisé au début du Modèle:S-. Les premières utilisations remontent à Huygens en 1654 et Takebe Kenkō en 1723, pour l'évaluation numérique de π.

Ce procédé est notamment utilisé pour définir une méthode numérique d'intégration : la méthode de Romberg, accélération de la méthode des trapèzes.

Présentation du principe

On suppose que la quantité inconnue Modèle:Mvar peut être approchée par une fonction Modèle:Math avec une convergence d'ordre Modèle:Mvar en Modèle:Mvar

$A-A(h)=a_n{h}^n+O(h^m),a_n\ne 0,m>n,$

expression dans laquelle le coefficient Modèle:Mvar n'est pas connu. Le principe d'extrapolation consiste à supprimer le terme en Modèle:Mvar par combinaison linéaire de deux valeurs de Modèle:Math, calculés avec des Modèle:Mvar différents : par exemple Modèle:Math et Modèle:Math. On obtient :

$R(h)=A(h/2)+\frac{A(h/2)-A(h)}{2^n-1}=\frac{2^{n}A(h/2)-A(h)}{2^n-1}$

Modèle:Math est l'extrapolé de Richardson qui approche Modèle:Mvar à l'ordre m>n en h. Le facteur 2 peut être remplacé par n'importe quel autre facteur. L'intérêt de la méthode est qu'il sera fréquemment plus aisé d'obtenir une précision donnée en calculant Modèle:Math que Modèle:Math avec un Modèle:Math beaucoup plus petit (risque d'erreur d'arrondi, grande quantité de calcul …).

Formule générale

On suppose que l'on dispose d'une approximation de Modèle:Mvar avec une formule d'erreur de cette forme

$A=A(h)+a_0{h}^{k_0}+a_1{h}^{k_1}+a_2{h}^{k_2}+\cdots a_z{h}^{k_z}+O(h^{k_{z+1}}),$

les coefficients étant inconnus. On se fixe un paramètre réel r>1 et on forme une combinaison entre la relation précédente et cette même relation prise au point $h/r$

$(r^{k_0}-1)A=r^{k_0}A\left(\frac{h}{r}\right)-A(h)+0+a_1(\frac{r^{k_0}}{r^{k_1}}-1)h^{k_1}+\dots +a_z(\frac{r^{k_0}}{r^{k_z}}-1)h^{k_z}+O(h^{k_{z+1}}).$

Le terme en Modèle:Math a disparu. Cette formule peut être itérée pour augmenter l'ordre, en évaluant Modèle:Math successivement aux points $h_n=\frac{h_0}{r^n}$.

Pour les formules d'erreur pouvant être exprimé sous la forme

$A=A(h)+a_1\varphi (h)+a_2{\varphi (h)}^2+a_3{\varphi (h)}^3+\cdots a_z{\varphi (h)}^z+O({\varphi (h)}^{z+1}),$

avec une fonction $\mathrm{\Phi }(x)$ connue telle que $\underset{x\to 0}{\lim }\varphi (x)=0$, un algorithme d'interpolation polynomial (par exemple l'algorithme d'Aitken-Neville) peut être utilisé. Dans ce cas, la suite des subdivisions Modèle:Mvar n'a pas nécessité d'être en progression géométrique.

Dans ce cas, le résultat de l'extrapolation de Richardson s'obtient en calculant la valeur en zéro du polynôme d'interpolation passant par les points $(\mathrm{\Phi }(h_i),A(h_i))$, où les Modèle:Mvar forment une suite décroissant vers 0. Une variante consiste à utiliser une interpolation par fraction rationnelle au lieu d'une interpolation polynomiale (Modèle:Lien ou le ρ-algorithme de P. Wynn).

Pour le cas encore plus général, où la formule d'erreur est de la forme :

$A=A(h)+a_1{\varphi }_1(h)+a_2{\varphi }_2(h)+a_3{\varphi }_3(h)+\cdots$

avec les fonctions ${\mathrm{\Phi }}_n(x)$ connues et telles que $\underset{x\to 0}{\lim }{\varphi }_n(x)=0$, l'algorithme d'interpolation linéaire généralisé de Mühlbach peut être utilisé : On obtient dans ce cas un algorithme d'accélération de la convergence appelé E-algorithme, de Håvie et Brezinski (à ne pas confondre avec l'ε-algorithme de P. Wynn).

Algorithme

L'algorithme présenté est celui associé à la méthode d'Aitken Neville, lorsque la formule d'erreur est de la forme :

$A=A(h)+a_1\varphi (h)+a_2{\varphi (h)}^2+a_3{\varphi (h)}^3+\cdots a_z{\varphi (h)}^z+O({\varphi (h)}^{z+1}),$

Pour N valeurs différentes calculés de Modèle:Mvar, avec la suite auxiliaire Modèle:Mvar correspondante, on initialise le tableau suivant :

$R_0^{(n)}=A_{n}n=0,1,\dots ,N-1$

Puis le reste du tableau (seulement une partie triangulaire) est calculé par les relations :

$R_{k+1}^{(n)}=R_k^{(n+1)}+\frac{R_k^{(n+1)}-R_k^{(n)}}{\frac{\varphi (h_n)}{\varphi (h_{n+k+1})}-1}\text{pour}n<N-1-k$

Le résultat de l'extrapolation de Richardson est $R_{N-1}^{(0)}$ (la pointe du triangle)

Le choix de la suite auxiliaire Modèle:Mvar est un compromis entre une suite décroissant suffisamment rapidement pour obtenir un algorithme numériquement stable, et une suite décroissant lentement évitant les difficultés de calcul des Modèle:Mvar. Pour des raisons de risque de divergence, on évitera d'utiliser une suite auxiliaire décroissant plus lentement que la suite harmonique (critère de Laurent).

D'autres algorithmes d'interpolation polynomiales peuvent être utilisés à la place de la méthode d'Aitken Neville, notamment la méthode de Lagrange (en particulier sa variante dite barycentrique). Dans ce cas précis, les polynômes formant la base de Lagrange peuvent être calculés une fois pour toutes (ceux-ci ne dépendant que de la suite Modèle:Mvar choisie), ce qui rend l'algorithme compétitif. Ceci est d'autant plus vrai s'il y a à extrapoler à plusieurs reprises, par exemple pour la résolution des équations différentielles, ou l'extrapolation d'une suite vectorielle.

Les suites classiquement utilisés sont :

• la suite de Romberg $h_n=\frac{h_0}{2^n}$
• les suites de Bulirsch $h_1=h_0/2;h_2=h_0/3;h_{n+2}=\frac{h_n}{2}$ ou bien $h_1=h_0/2;h_2=h_0/4;h_3=h_0/6;h_{n+2}=\frac{h_n}{2}$
• la suite de Deuflhard $h_n=\frac{h_0}{2n}$

Lorsque les exposants de la formule d'erreur de Modèle:Math ne suivent pas une progression arithmétique, les algorithmes précédents doivent être adaptés, avec la contrainte d'utiliser une suite Modèle:Mvar géométrique, du type Modèle:Math (par exemple, celle de Romberg, de raison 2). D'autre part, les différents exposants de la formule d'erreur doivent être connus, car ils rentre dans la formule de calcul du tableau.

Exemples d'application

Approximation numérique de la dérivée d'une fonction

Le point de départ est le développement de Taylor de la fonction à dériver :

$f(x+h)=f(x)+f^{\prime }(x)h+\frac{f^{″}(x)}{2!}{h}^2+\frac{f^{‴}(x)}{3!}{h}^3+\frac{f^{⁗}(x)}{4!}{h}^4+\cdots$

la dérivée de Modèle:Math est déduite

$f^{\prime }(x)=\frac{f(x+h)-f(x)}{h}-\frac{f^{″}(x)}{2!}h-\frac{f^{‴}(x)}{3!}{h}^2-\frac{f^{⁗}(x)}{4!}{h}^3+\cdots .$

en appelant $A(h)=\frac{f(x+h)-f(x)}{h}$, le terme d'erreur se présente sous une forme :

$A(h)=f^{\prime }(x)+a_{1}h+a_2{h}^2+a_3{h}^3+\cdots$

forme convenable pour être accélérée à l'aide de l'extrapolation de Richardson. Il suffira de calculer Modèle:Math pour différentes valeur de Modèle:Mvar (par exemple la suite de Romberg), et de s'en servir pour éliminer les termes en Modèle:Mvar, Modèle:MathModèle:Etc.

En soustrayant le développement de Taylor de Modèle:Math à celui de Modèle:Math, on obtient :

$f(x+h)-f(x-h)=2f^{\prime }(x)h+2\frac{f^{‴}(x)}{3!}{h}^3+2\frac{f^{(5)}(x)}{5!}{h}^5+\cdots$

on en tire une autre approximation de la dérivée :

$A(h)=\frac{f(x+h)-f(x-h)}{2h}=f^{\prime }(x)+a_2{h}^2+a_4{h}^4+a_6{h}^6+\cdots .$

qui ne contient que des termes d'erreur de degré pair. L'utilisation de l'extrapolation de Richardson est dans ce cas beaucoup plus efficace, car à chaque étape, le degré d'approximation est augmenté de 2.

De la même manière, cette fois-ci en additionnant les développements de Modèle:Math et Modèle:Math, on en tire une formule du même type servant à calculer la dérivée seconde :

$A(h)=\frac{f(x+h)-2f(x)+f(x-h)}{h^2}=f^{″}(x)+a_2{h}^2+a_4{h}^4+a_6{h}^6+\cdots .$

Dans les deux cas précédents, on pourra utiliser la méthode de Aitken-Neville pour le calcul de l'extrapolation de Richardson, avec la fonction $\mathrm{\Phi }(h)=h^2$.

Intégration numérique

La méthode des trapèzes possède un terme d'erreur n'ayant que des puissances paires de h, d'après la formule d'Euler-Maclaurin. Elle est donc très bien adaptée à l'utilisation de l'extrapolation de Richardson. La méthode résultante est appelée méthode de Romberg.

La méthode du point médian possède aussi ce type de développement, et peut donc fournir une variante de la méthode de Romberg. Les deux méthodes peuvent être utilisées de concert pour fournir une estimation de l'erreur de l'intégrale résultante.

Intégration numérique des équations différentielles

La méthode d'intégration de base qui se prête le mieux à l'extrapolation de Richardson est la méthode du point milieu (où du point milieu modifié) de Gragg. Elle possède aussi un terme d'erreur où n'apparait que les puissances paires de Modèle:Mvar (lorsque les subdivisions sont paires). La méthode résultante est connue sous le nom de Modèle:Lien. Il existe plusieurs variantes, dépendant du choix de la séquence des Modèle:Mvar, du type d'extrapolation (polynomiale ou rationnelle), et du calcul de la taille du pas global.

Voir aussi

Article connexe

Delta-2