Epsilon algorithme

Pour les articles homonymes, voir Epsilon (homonymie).

En analyse numérique, l'ε-algorithme est un algorithme non linéaire d'accélération de la convergence de suites numériques. Cette algorithme a été proposé par Peter Wynn^[1] en 1956 pour calculer la transformation de Shanks. C'est l'un des algorithmes d'accélération de la convergence les mieux connus et des plus utilisés. C'est une généralisation de l'algorithme Delta-2 d'Aitken.

Présentation

Soit une suite numérique S_n dont on cherche à connaitre la limite S. l'ε-algorithme consiste à calculer un tableau en initialisant la première colonne par la suite S_n, et en calculant les autres cases à l'aide des relations suivantes :

$\epsilon^{(n)}_{-1}=0 \text{ }\epsilon^{(-n-1)}_{2n}=0 \text{ } \epsilon^{(n)}_{0}=S_n \text{ }\forall n$

$\epsilon^{(n)}_{k+1}= \epsilon^{(n+1)}_{k-1} + \frac{1}{\epsilon^{(n+1)}_{k} - \epsilon^{(n)}_{k}} \text{ }\forall n,k$

Les différentes colonnes d'indice k paires du tableau fournissent d'autres suites convergeant en général plus vite que la suite S_n d'origine. Les colonnes d'indice impaires peuvent êtres considérés comme des intermédiaires de calcul.

On présente traditionnellement les résultats de l'ε-algorithme sous forme d'un tableau aux lignes décalés: la formule de calcul de l'ε-algorithme relie ainsi les termes formant un losange dans le tableau.

$\begin{matrix} 0 & & 0 & & 0 & \cdots \\ & \epsilon^{(0)}_{0}=S_0 & & \epsilon^{(-1)}_{2} & & \epsilon^{(-2)}_{4} \cdots \\ 0 & & \epsilon^{(0)}_{1} & & \epsilon^{(-1)}_{3} & \cdots \\ & \epsilon^{(1)}_{0}=S_1 & & \epsilon^{(0)}_{2} & & \epsilon^{(-1)}_{4} \cdots \\ 0 & & \epsilon^{(1)}_{1} & & \epsilon^{(0)}_{3} & \cdots \\ & \epsilon^{(2)}_{0}=S_2 & & \epsilon^{(1)}_{2} & & \epsilon^{(0)}_{4} \cdots \\ 0 & & \epsilon^{(2)}_{1} & & \epsilon^{(1)}_{3} & \cdots \\ & \epsilon^{(3)}_{0}=S_3 & & \epsilon^{(2)}_{2} & & \epsilon^{(1)}_{4} \cdots \\ & \vdots & & \vdots & & \vdots \\ \end{matrix}$

Propriétés

Le tableau obtenu par l'ε-algorithme est une des méthodes utilisées pour calculer la transformation de Shanks d'une suite. En effet:

$e_k(S_n)= \epsilon^{(n)}_{2k} \text{ }\forall k,n.$

où e_k(S_n) est la transformée de Shanks d'ordre k de la suite S_n. L'ε-algorithme partage donc toutes les propriétés de la transformée de Shanks. Notamment, l'ε-algorithme est un moyen pour calculer la table de Padé d'une série entière. De même, la suite ε⁽ⁿ⁾₂ correspond au Delta-2 d'Aitken. Les approximants de Padé à plusieurs variables peuvent aussi très simplement être calculés à l'aide de l'ε-algorithme. Par exemple, pour une fonction à 2 variables, dont le développement de Taylor est:

$f(x,y)= \sum_{i=0}^{\infty} \sum_{j=0}^{\infty} a_{i,j}.x^i.y^j$

Sa table de Padé peut être obtenue à l'aide de l'ε-algorithme, en initialisant la première colonne par:

$\epsilon^{(n)}_{0}=\sum_{m=0}^{n}\left(\sum_{i+j=m} a_{i,j}.x^i.y^j \right)$

Variantes

ε-algorithme vectoriel

Il existe plusieurs variantes de l'ε-algorithme pouvant êtres utiliser avec des suites vectorielles (ε-algorithme vectoriel, topologique ou scalaire appliqué à chacune des composantes du vecteur). L'ε-algorithme vectoriel est analogue dans sa forme à l'ε-algorithme classique:

$\epsilon^{(n)}_{-1}=0 \text{ }\epsilon^{(-n-1)}_{2n}=0 \text{ } \epsilon^{(n)}_{0}=S_n \text{ }\forall n$

$\epsilon^{(n)}_{k+1}= \epsilon^{(n+1)}_{k-1} + {(\epsilon^{(n+1)}_{k} - \epsilon^{(n)}_{k})}^{-1} \text{ }\forall n,k$

où les suites S_n et les 0 sont des vecteurs.

Reste à définir ce qu'est l'inverse d'un vecteur! On pourra prendre par exemple comme inverse d'un vecteur U:

$U^{-1}=\frac {U}{{\|U\|}^2}$

où la norme du vecteur est la norme euclidienne. L'ε-algorithme vectoriel a de nombreuses applications pratiques mais est très sensible aux erreurs d'arrondi. Il permet notamment de généraliser la méthode de Aitken-Steffensen à des systèmes d'équations non linéaires, fournissant ainsi un algorithme à convergence quadratique, ne nécessitant pas de calcul du Jacobien, à l'inverse de la méthode de Newton.

L'ε-algorithme matriciel existe pour les matrices carrés, dans ce cas l'inverse apparaissant dans l'algorithme est tout simplement celui d'une matrice carré classique.

ε-algorithme confluent

L'ε-algorithme classique accélère une suite numérique, c'est-à-dire une fonction d'une variable discrète 'n'. La version confluente de cet algorithme accélère la convergence d'une fonction d'une variable continue 't'. On l'obtient en faisant un changement de variable x=t+n*h à l'algorithme d'origine, et en faisant tendre h vers 0. on obtient:

$\epsilon_{-1}(t)=0 \text{ } \epsilon_{0}(t)=f(t)$

$\epsilon_{k+1}(t)= \epsilon_{k-1}(t) + \frac{1}{\epsilon'_{k}(t)} \text{ } k=0,1,...$

Cet algorithme est utilisé par exemple pour l'évaluation des intégrales impropres lorsque les dérivées de la fonction sont accessibles, ou à l'élaboration de séries asymptotiques pour certaines fonctions.

La dérivée des termes apparaissant dans la formule peuvent se ramener aux dérivées de la fonction de base f(t) par calcul formel, où en utilisant les relations des déterminants de Henkel:

$\epsilon'_{2k}(t)=\frac{H_{k}^{(1)} H_{k+1}^{(1)}}{(H_{k}^{(2)})^2}$

$\epsilon'_{2k+1}(t)= -\frac{H_{k}^{(2)} H_{k+1}^{(2)}}{(H_{k+1}^{(1)})^2}$

avec

$H_{k+2}^{(n-1)} H_{k}^{(n+1)}+(H_{k+1}^{(n)})^2 = H_{k+1}^{(n-1)} H_{k+1}^{(n+1)}\text{ } n,k=0,1,...$

et

$H_{0}^{(n)}=1 \text{ } H_{1}^{(n)}=f^{(n)}(t)\text{ } n=0,1,...$

On trouve, en explicitant les premières expressions:

$\frac{1}{\epsilon'_{1}(t)}=- \frac{[f'(t)]^2}{f''(t)}$

$\frac{1}{\epsilon'_{2}(t)}=- \frac{f'(t)f'''(t)-f''^2(t)}{f''(t)[f''(t)f^{(4)}(t)-f'''^2(t))]^2}$

ε-algorithme dépendant d'une suite auxiliaire

Il arrive fréquemment que la suite à accélérer S_n dépende d'une suite auxiliaire x_n, celle-ci tendant vers l'infini (par exemple des approximations par discrétisation avec un maillage de plus en plus fin). Il est possible d'utiliser l'ε-algorithme de base dans ces cas mais la façon dont évolue la suite x_n associée (par exemple l'évolution du nombre de mailles entre chaque S_n) est un renseignement précieux qui pourrait aider l'accélération de la convergence. Deux variantes de l'ε-algorithme utilise cet information supplémentaire et donne souvent de meilleurs résultats.

$\epsilon^{(n)}_{-1}=0 \text{ }\epsilon^{(-n-1)}_{2n}=0 \text{ } \epsilon^{(n)}_{0}=S_n \text{ }\forall n$

$\epsilon^{(n)}_{k+1}= \epsilon^{(n+1)}_{k-1} + \frac{x_{n+1}-x_n}{\epsilon^{(n+1)}_{k} - \epsilon^{(n)}_{k}} \text{ }\forall n,k$

Et la deuxième variante:

$\epsilon^{(n)}_{-1}=0 \text{ }\epsilon^{(-n-1)}_{2n}=0 \text{ } \epsilon^{(n)}_{0}=S_n \text{ }\forall n$

$\epsilon^{(n)}_{k+1}= \epsilon^{(n+1)}_{k-1} + \frac{x_{n+k+1}-x_{n+k}}{\epsilon^{(n+1)}_{k} - \epsilon^{(n)}_{k}} \text{ }\forall n,k$

Ces algorithmes s'apparentent aux méthodes par extrapolation (extrapolation de Richardson ou ρ-algorithme). Ils peuvent accélérer des suites particulièrement récalcitrantes avec d'autres algorithmes.

ε-algorithme d'interpolation

L'ε-algorithme classique permet de calculer la table de Padé d'une fonction lorsque les termes de la suite S_n sont les sommes partielles du développement limité de cette fonction. Par analogie, il est possible aussi de construire la table des fractions rationnelles d'interpolation d'une fonction en partant de la suite des polynômes d'interpolation.

$\epsilon^{(n)}_{-1}=0 \text{ }\epsilon^{(-n-1)}_{2n}=0 \text{ } \epsilon^{(n)}_{0}=P_n(x) \text{ }\forall n$

$\epsilon^{(n)}_{k+1}= \epsilon^{(n+1)}_{k-1} + \frac{1}{(x-x_{n+k+1})\cdot (\epsilon^{(n+1)}_{k} - \epsilon^{(n)}_{k})} \text{ }\forall n,k$

$R_{n,k}(x)= \epsilon^{(n-k)}_{2k} \text{ }\forall k,n.$

P_n(x) étant le polynôme d'interpolation passant par les points (x₀,y₀),(x₁,y₁),...,(x_n,y_n), que l'on pourra calculer par la méthode de Lagrange, de de Newton où de Neville, par exemple. R_n,k(x)est la fraction rationnelle de degré n au numérateur et k au dénominateur passant par les points (x₀,y₀),(x₁,y₁),...,(x_n+k,y_n+k)

Cet algorithme peut se présenter comme un complément aux méthodes de Thiele ou Bulirsch-Stoer, avec la particularité que l'ε-algorithme permet de calculer tout le tableau des fractions rationnelles d'interpolation au lieu d'une seule diagonale, ou un escalier. De plus, si le polynôme utilisé est un polynôme d'interpolation d'Hermite (imposant d'interpoler la fonction, mais aussi sa ou ses dérivées), l'ε-algorithme d'interpolation fournira les fraction rationnelle d'interpolation vérifiant les mêmes critères que le polynôme d'Hermite. La méthode d'interpolation de Newton avec confluence des abscisses (pour les points où les dérivées sont à renseigner) est toute indiquée pour calculer le polynôme d'Hermite, d'autant que la suite des abscisses utilisée (avec les confluences) sera nécessaire pour l'ε-algorithme d'interpolation.

Règles particulières

Il peut arriver lors du calcul du tableau de l'ε-algorithme, qu'à certains endroits se produisent des divisions par 0 où presque. Ceci peut générer des dégradations de précision pour le calcul des zones du tableau dépendant de la case suspecte, voir des plantage de programme. Plusieurs auteurs ont développés des algorithmes spéciaux^[2] pour contourner les cases problématiques en cas de détection de risque de division par 0. Ces algorithmes ne fonctionnent que si les cases singulières ne sont pas trop proches les unes des autres.

Exemples

Voir aussi les exemples de la transformée de Shanks.

• accélération de la convergence d'une suite

La méthode de Bernoulli permet, sous certaines conditions, d'évaluer la plus grande (ou la plus petite) racine d'un polynôme donné. Mais la suite générée par la méthode peut converger lentement vers la valeur de la racine: l'ε-algorithme permet d'accélérer cette convergence. Par exemple, avec le polynôme x²-x-1 dont les racines sont le nombre d'or φ et 1/φ, la suite générée par la méthode de Bernoulli donne la suite de Fibonacci F_n dont le ratio des termes successifs converge vers φ=1,6180339887499...

$\scriptstyle{\begin{matrix} & \epsilon^{(n)}_{0}=F_{n+1}/F_n & & \epsilon^{(n)}_{2} & & \epsilon^{(n)}_{4} & & \epsilon^{(n)}_{6} & & \epsilon^{(n)}_{8 } \\ 0 \\ & 1 \\ 0 & & -2 \\ & 2 & & 1,62500000 \\ 0 & & 6 & & -162 \\ & 1,5 & & 1,61904761 & & 1,61805555 \\ 0 & & -15 & & -1170 & & -45090,00 \\ & 1,66666667 & & 1,61818181 & & 1,61803278 & & 1,6180339985 \\ 0 & & 40 & & -7880 & & 780240,00 & & -103779601,9 \\ & 1,60000000 & & 1,61805555 & & 1,61803405 & & 1,6180339889 & & 1,61803398875 \\ 0 & & -104 & & -54392 & & -14196624,02 & & -4926280210,4 \\ & 1,62500000 & & 1,61803713 & & 1,61803398 & & 1,6180339888 \\ 0 & & 273 & & -371826 & & 253400249,03 \\ & 1,61538461 & & 1,61803444 & & 1,61803398 \\ 0 & & -714 & & -2551122 \\ & 1,61904761 & & 1,61803405 \\ 0 & & 1870 \\ & 1,61764705 \\ \end{matrix}}$

Dans cet exemple, seuls les termes d'indice haut positif ont été calculés.

Nous constatons dans cette exemple qu'effectivement seules les colonnes paires convergent vers la suite d'origine, et ceci plus rapidement (12 chiffres exacts au lieu de 3 dans la suite initiale).

Il est à noter que l'ε-algorithme possède une propriété particulière vis-à-vis de la méthode de Bernoulli qui permet de l'utiliser aussi dans un autre but que l'accélération de la convergence: le calcul simultané de toutes les racines du polynôme au lieu d'une seule à la fois.

Soit le polynôme:

P(x) = a₀ + a₁.x + a₂.x² + ... + a_p.x^p

dont les racines λ_k sont réelles, telles que |λ₁|> |λ₂|>...>|λ_p|

la méthode de Bernoulli génère la suite y_n définie par:

$y_{n+p}= \frac{-1}{a_p}\left( a_{p-1} . y_{n+p-1}+...+a_1 . y_{n+1} + a_0 . y_{n} \right)\text{ }n=0,1...$

dont le rapport des termes consécutifs y_n+1/y_n converge vers λ₁.

La suite est initialisée avec y₀, y₁...y_p-1 arbitraires, non tous nuls.

en appliquant l'ε-algorithme , non pas sur le ratio y_n+1/y_n, mais directement sur la suite des y_n, on a:

$\lim_{n \to \infty} \frac{\epsilon^{(n+1)}_{2k}}{\epsilon^{(n)}_{2k}}=\lambda_{k+1} \text{ }k=0,1...p-1$

et

$\lim_{n \to \infty} \frac{\epsilon^{(n)}_{2k+1}}{\epsilon^{(n+1)}_{2k+1}}=\lambda_{k+1} \text{ }k=0,1...p-1$

On constate dans ce cas que les colonnes impaires s'avèrent aussi intéressantes.

la vitesse de convergence, pour k fixé est:

$\frac{\epsilon^{(n+1)}_{2k}}{\epsilon^{(n)}_{2k}}=\lambda_{k+1}+ O[(\frac{\lambda_{k+2}}{\lambda_{k+1}})^{n+k+1}]$

On peut donc utiliser à nouveau l'ε-algorithme, mais cette fois pour accélérer la convergence de chacun de ces ratio, comme on l'a fait pour la suite y_n+1/y_n.

• Résolution d'un système d'équation non linaire

L'une des utilisation de l'ε-algorithme est de fournir une généralisation de la méthode d'Aitken-Steffensen pour des systèmes d'équations non linéaires.

La méthode consiste à réécrire le système d'équation F(X)=0 en un système de la forme X=G(X) (X, F et G étant des vecteurs représentant les inconnues ou le système d'équations). En partant d'un vecteur X₀ arbitraire (de préférence proche de la solution du système) générer la suite de vecteurs X_k de la manière suivante:

- U₀ = X_k k = 0,1...

- U_{n + 1} = G(U_n) n = 0,1...2p

- calculer l'ε-algorithme de cette suite de vecteurs U_n

- $X_{k+1}=\epsilon^{(0)}_{2p}$

p étant le nombre d'équations du système.

La suite des X_k générés par cette méthode est à convergence quadratique.

Plusieurs choix pour G étant possibles, il est préférable, bien que non nécessaire, de choisir G de telle sorte que la suite généré U_n+1=G(U_n) converge vers la solution du système.

L'ε-algorithme vectoriel parait le plus adapté à cet algorithme, mais l'ε-algorithme scalaire appliqué à chacune des composantes du vecteur fonctionne aussi.

• Fraction rationnelle d'interpolation

Cette section est vide, insuffisamment détaillée ou incomplète. Votre aide est la bienvenue !

• Variante des méthodes utilisant l'extrapolation de Richardson

Plusieurs méthodes numériques utilisent l'extrapolation de Richardson pour accélérer la convergence de méthodes d'ordre peu élevé. On trouve par exemple son utilisation dans:

- l'évaluation de la dérivée d'une fonction, par accélération de la convergence de la méthode des différences centrales

- l'intégration numérique, en accélérant la méthode des trapèzes (méthode de Romberg) ou la méthode du point médian

- l'intégration des équation différentielles, en accélérant la méthode du point milieu (méthode de Gragg-Bulirsch-Stoer)

Dans chacune de ces applications, il est possible de remplacer l'extrapolation de Richardson par l'ε-algorithme ou ses généralisations. La suite associée x_n devant tendre vers l'infini dans le cas des généralisation de l'ε-algorithme, on prendra l'inverse de la suite associée à l'extrapolation de Richardson. Cependant l'extrapolation de Richardson étant particulièrement bien adaptée à accélérer les suites générées par ces méthodes, l'ε-algorithme en général est moins performant. Dans certains cas, par exemple le calcul numérique d'une intégrale impropre, l'extrapolation de Richardson échoue à accélérer la convergence, et l'ε-algorithme devient alors la méthode de choix.

Notes

1. ↑ P. Wynn, On a device for computing the e_m(S_n) transformation, MTAC, 10(1956), p91-96
2. ↑ Wynn Peter, « Singular rules for certain non-linear algorithms », dans BIT, vol. 3, 1963, p. 175–195 [texte intégral, lien DOI] 

Références

• Claude Brezinski, Algorithmes d'accélération de la convergence: étude numérique, éditions Technip, 1978, (ISBN 2-7108-0341-0)

Voir aussi

Articles connexes

• Transformation de Shanks
• Delta-2
• ρ-algorithme
• Extrapolation de Richardson