Extrapolation de Richardson

On suppose que la quantité inconnue A peut être approchée par une fonction A(h) avec une convergence d'ordre n en h

${\displaystyle A-A(h)=a_{n}h^{n}+O(h^{m}),~a_{n}\neq 0,~m>n,}$

expression dans laquelle le coefficient a_n n'est pas connu. Le principe d'extrapolation consiste à supprimer le terme en hⁿ par combinaison linéaire de deux valeurs de A(h), calculés avec des h différents : par exemple A(h) et A(h/2). On obtient :

${\displaystyle 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}}}$

R(h) est l'extrapolé de Richardson qui approche A à 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 R(h) que A(h') avec un h' beaucoup plus petit (risque d'erreur d'arrondi, grande quantité de calcul …).

On suppose que l'on dispose d'une approximation de A avec une formule d'erreur de cette forme

${\displaystyle 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 ${\displaystyle h/r}$

${\displaystyle {(r^{k_{0}}-1)A=r^{k_{0}}A\left({\frac {h}{r}}\right)-A(h)+0+a_{1}\left({\frac {r^{k_{0}}}{r^{k_{1}}}}-1\right)h^{k_{1}}+\dots +a_{z}\left({\frac {r^{k_{0}}}{r^{k_{z}}}}-1\right)h^{k_{z}}+O(h^{k_{z+1}}).}}$

Le terme en h^k₀ a disparu. Cette formule peut être itérée pour augmenter l'ordre, en évaluant A(h) successivement aux points ${\displaystyle h_{n}={\frac {h_{0}}{r^{n}}}}$.

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

${\displaystyle A=A(h)+a_{1}{\phi (h)}+a_{2}{\phi (h)}^{2}+a_{3}{\phi (h)}^{3}+\cdots a_{z}{\phi (h)}^{z}+O({\phi (h)}^{z+1}),}$

avec une fonction ${\displaystyle \Phi (x)}$ connue telle que ${\displaystyle \lim _{x\to 0}\phi (x)=0}$, un algorithme d'interpolation polynomial (par exemple l'algorithme d'Aitken-Neville) peut être utilisé. Dans ce cas, la suite des subdivisions h_n 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 ${\displaystyle (\Phi (h_{i}),A(h_{i}))}$, où les h_i forment une suite décroissant vers 0. Une variante consiste à utiliser une interpolation par fraction rationnelle au lieu d'une interpolation polynomiale (Bulirsch et Stoer (en)[4] ou le ρ-algorithme de P. Wynn).

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

${\displaystyle A=A(h)+a_{1}\phi _{1}(h)+a_{2}\phi _{2}(h)+a_{3}\phi _{3}(h)+\cdots }$

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

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

${\displaystyle A=A(h)+a_{1}{\phi (h)}+a_{2}{\phi (h)}^{2}+a_{3}{\phi (h)}^{3}+\cdots a_{z}{\phi (h)}^{z}+O({\phi (h)}^{z+1}),}$

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

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

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

${\displaystyle R_{k+1}^{(n)}=R_{k}^{(n+1)}+{\frac {R_{k}^{(n+1)}-R_{k}^{(n)}}{{\frac {\phi (h_{n})}{\phi (h_{n+k+1})}}-1}}\quad {\textrm {pour}}\quad n<N-1-k}$

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

Le choix de la suite auxiliaire h_n 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 A_n. 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 h_n 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 ${\displaystyle h_{n}={\frac {h_{0}}{2^{n}}}}$
• les suites de Bulirsch ${\displaystyle h_{1}=h_{0}/2;h_{2}=h_{0}/3;h_{n+2}={\frac {h_{n}}{2}}}$ ou bien ${\displaystyle 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 ${\displaystyle h_{n}={\frac {h_{0}}{2n}}}$

Lorsque les exposants de la formule d'erreur de A(h) ne suivent pas une progression arithmétique, les algorithmes précédents doivent être adaptés, avec la contrainte d'utiliser une suite h_n géométrique, du type h_n=h₀/rⁿ (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.