Epsilon algorithme

En analyse numérique, l'ε-algorithme est un algorithme non linéaire d'accélération de la convergence de suites numériques. Cet 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 les plus utilisés. C'est une généralisation de l'algorithme Delta-2 d'Aitken.

Présentation[modifier | modifier le code]

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

Les différentes colonnes d'indice k pair du tableau fournissent d'autres suites convergeant en général plus vite que la suite Sn d'origine. Les colonnes d'indice impair peuvent être considérées 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ées : la formule de calcul de l'ε-algorithme relie ainsi les termes formant un losange dans le tableau.

Propriétés[modifier | modifier le code]

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

ek(Sn) est la transformée de Shanks d'ordre k de la suite Sn. 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 ε(n)2 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 à deux variables, dont le développement de Taylor est :

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

Variantes[modifier | modifier le code]

ε-algorithme vectoriel[modifier | modifier le code]

Il existe plusieurs variantes de l'ε-algorithme pouvant être utilisées 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 :

où les suites Sn 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 :

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. Il 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ées, dans ce cas l'inverse apparaissant dans l'algorithme est tout simplement celui d'une matrice carrée classique.

ε-algorithme confluent[modifier | modifier le code]

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+nh à l'algorithme d'origine, et en faisant tendre h vers 0. on obtient :

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 :

avec

et

On trouve, en explicitant les premières expressions :

ε-algorithme dépendant d'une suite auxiliaire[modifier | modifier le code]

Il arrive fréquemment que la suite à accélérer Sn dépende d'une suite auxiliaire xn, 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 xn associée (par exemple l'évolution du nombre de mailles entre chaque Sn) est un renseignement précieux qui pourrait aider l'accélération de la convergence. Deux variantes de l'ε-algorithme utilisent cette information supplémentaire et donnent souvent de meilleurs résultats.

Et la deuxième variante :

Ces algorithmes s'apparentent aux méthodes par extrapolation (extrapolation de Richardson ou ρ-algorithme). Ils peuvent accélérer des suites à convergence logarithmique, que l'ε-algorithme classique n'accélère pas.

ε-algorithme d'interpolation[modifier | modifier le code]

L'ε-algorithme classique permet de calculer la table de Padé d'une fonction lorsque les termes de la suite Sn sont les sommes partielles du développement limité de cette fonction. Une variante[2] permet de construire la table des fractions rationnelles d'interpolation d'une fonction en partant de la suite des polynômes d'interpolation.

Le polynôme Pn(x) étant la valeur en x du polynôme d'interpolation passant par les points (x0,y0) , (x1,y1),...,(xn,yn), que l'on pourra calculer par la méthode de Lagrange, de Newton où de Neville, par exemple. Rn,k(x) est la fraction rationnelle de degré n au numérateur et k au dénominateur passant par les points (x0,y0) , (x1,y1),...,(xn+k,yn+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[modifier | modifier le code]

Il peut arriver lors du calcul du tableau de l'ε-algorithme, qu'à certains endroits se produisent des divisions par 0 ou 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 plantages de programme. Plusieurs auteurs ont développé des algorithmes spéciaux[3] 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[modifier | modifier le code]

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

Accélération de la convergence d'une suite[modifier | modifier le code]

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 x2x – 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 Fn dont le ratio des termes successifs converge vers φ=1,6180339887499...

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

Nous constatons dans cet 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 :

dont les racines λk sont réelles, telles que |λ1|> |λ2|>...>|λp|

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

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

La suite est initialisée avec y0, y1,... yp–1 arbitraires, non tous nuls.

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

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

La vitesse de convergence, pour k fixé, est :

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

Résolution d'un système d'équation non linéaire[modifier | modifier le code]

L'une des utilisations 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 X0 arbitraire (de préférence proche de la solution du système), on génère la suite de vecteurs Xk de la manière suivante :

  • calculer l'ε-algorithme de cette suite de vecteurs Un

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

La suite des Xk 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ée Un+1=G(Un) 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[modifier | modifier le code]

Exemple: interpolation polynomiale de la fonction logarithme népérien et calcul des fraction rationnelles correspondantes à l'aide de l'ε-algorithme d'interpolation

La fonction logarithme sera interpolée en trois points: 1/2, 1 et e. Pour tester la mise en œuvre de l'algorithme dans le cas d'une interpolation d'Hermite, l'interpolation en 1/2 sera double (interpolation du logarithme et de sa dérivée), et triple en 1 (interpolation du logarithme, et de ses dérivées premières et secondes). Il s'agit d'un polynôme de degré 5 qui doit donc vérifier :

Tableau des différences divisées pour l'interpolation polynomiale de ln(x) en 3 points
i xi yi / d0 d1 d2 d3 d4 d5
0 2,71828182845904 1 0,581976706869326 -0,243279819530861 0,149405165216329 -0,178413885100506 0,24817470934455
1 1 0 1 -0,5 0,545177444479562 -0,728935333122626
2 1 0 1 -0,772588722239781 0,909645111040875
3 1 0 1,38629436111989 -1,22741127776022
4 0.5 -0,693147180559945 2
5 0.5 -0,693147180559945

Les cases en italiques correspondent à des indéterminées des différences divisées (du fait des points doubles ou triples), et sont remplacés par les valeurs de la dérivée en ce point (ou de f"(x)/2! pour la différence seconde). Les cases en gras dans le tableau sont celles utiles au calcul de la forme de Newton du polynôme d'interpolation. Le polynôme s'écrit, à l'aide de la formule de Newton :

.
  • Calcul de P5(2), et application de l'ε-algorithme d'interpolation.

La suite des Pi(2) servant à initialiser le tableau de l'ε-algorithme d'interpolation, est la suite de polynômes de degré 0, 1, 2... obtenue en retenant 1,2,3... termes du polynôme P5(2) écrit sous la forme de Newton ci dessus. Cette suite de polynôme passent respectivement par les points d'abscisse (x0,y0); (x0,y0),(x1,y1); (x0,y0),(x1,y1),(x2,y2); ...

Comparaison de l'interpolation polynomiale et rationnelle par l'ε-algorithme d'interpolation
Interpolation polynomiale de ln(2)=0,693147180559945..., et calcul de l'ε-algorithme d'interpolation
i xi ε(i)-1 Pi(x)= ε(i)0 ε(i)1 ε(i)2 ε(i)3 ε(i)4
0 2,71828182845904 0 1 -2,39221119117733 0,705207033512282 -61,0702755830021 0,693879569827545
1 1 0 0,581976706869326 5,72267438319407 0,69023539353372 121,870014050176 0,692833802968095
2 1 0 0,756720180469139 -9,31836050756016 0,695317144633344 -146,585461084686
3 1 0 0,649405165216329 5,20217803449192 0,690925043467957
4 0.5 0 0,777556616828803 -2,49324571003365
5 0.5 0 0,51016754082086


La valeur de P5(2), pourtant située dans l'intervalle d'interpolation, est assez éloignée de ln(2), valeur de la fonction à interpoler. Sur le graphique, on note une forte dégradation des qualités d'interpolation du polynôme entre 1.5 et e. La présence d'un pôle à proximité de l'intervalle d'interpolation, ainsi que la croissance logarithmique de la courbe se prête mal à l'interpolation par polynômes. Ceci se prête davantage à l'interpolation rationnelle. L'application de l'ε-algorithme d'interpolation à ce polynôme améliore significativement les qualités d'interpolation sur cet exemple. Sur le graphique, certaines fractions rationnelles obtenues sont presque indiscernable de la fonction à interpoler, dans l'intervalle d'interpolation, et au delà. Les nombres en gras dans la table de l'ε-algorithme représente les fraction rationnelles d'interpolation de degré 5/0, 4/1 et 3/2. Elle peut être continuée pour intégrer celles de degré 2/3, 1/4 et 0/5.

Accélération de la convergence des séries de Fourier et de Tchebychev[modifier | modifier le code]

Les sommes partielles des séries de Fourier ou de Tchebychev peuvent fournir une suite qui peut être accélérée par l'ε-algorithme. Cependant, une manière plus efficace de procéder est de transformer, par un changement de variable, ces séries en une série entière. Les sommes partielles de ces séries entières seront les données d'entrée pour l'ε-algorithme.

La série de Fourier :

ou celle de Tchebychev:

se transforment chacune en une série entière (d'une variable complexe)

ou

dont la partie réelle est la série d'origine. Les Tn(x) représentent les polynômes de Tchebychev de 1ère espèce.

Ces séries se prêtent idéalement à l'accélération de la convergence par l'ε-algorithme (avec des nombres complexes) sur leurs sommes partielles. La partie réelle du résultat sera la valeur accélérée de la série de Fourier ou de Tchebychev d'origine.

Par exemple, la fonction Logarithme sur [1, 9], dont la série de Tchebychev vaut[4] :

devient après changement de variable la série entière :

Le graphique ci contre compare la fonction Logarithme avec sa série de Tchebychev en retenant les termes jusqu'au degré 6, ainsi que la valeur accélérée de cette série partielle par l'ε-algorithme. On note une nette amélioration de la convergence (erreur divisée par 80 ou plus). Le résultat accéléré conserve les caractéristiques typiques des approximations de Tchebychev: les oscillations d'amplitudes approximativement égales de l'erreur sur l'intervalle concerné (à peine perceptible sur le graphique).

Accélération de la convergence de série de Tchebychev de ln(x) par l ε-algorithme

Formellement, l'ε-algorithme transforme la série auxiliaire en approximant de Padé. En explicitant les calculs, l'approximant de Padé ici vaut:

dont la partie réelle, en revenant à la variable d'origine, vaut:

Cette fraction a les caractéristiques proches de la fraction min-max, optimale pour l'approximation: elle peut fournir un bon point de départ pour initialiser l'algorithme de Remez rationnel.


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 équations 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 xn devant tendre vers l'infini dans le cas des généralisations 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[5].

Par exemple, l'intégrale

possède une singularité algébrique en zéro, qui ralentit la convergence de la méthode de Romberg. L'utilisation de l'ε-algorithme à la place de l'extrapolation de Richardson dans cette méthode donne de meilleurs résultats :

Méthode de Romberg et l'ε-algorithme avec une intégrale impropre
subdivisions Trapèzes Romberg ε-algorithme
1 0,5 0,5 0,5
2 0,603553390593274 0,638071187457699 0,603553390593274
4 0,643283046242747 0,657756603281563 0,668014371343285
8 0,658130221624454 0,663607569112292 0,666989411481855
13 0,663581196877228 0,665592865129466 0,666661304912164
32 0,665558936278942 0,666287699033842 0,666666280204189
64 0,666270811378507 0,666532741199894 0,666666671581115

Notes et références[modifier | modifier le code]

  1. (en) P. Wynn, « On a device for computing the em(Sn) transformation », MTAC, 10 (1956), p. 91-96.
  2. (en) Guido Claessens, Some aspects of the rational Hermite interpolation table and its applications : Ph. D. thesis, Université d'Anvers, .
  3. (en) Peter Wynn, « Singular rules for certain non-linear algorithms », BIT, vol. 3,‎ , p. 175–195 (DOI 10.1007/BF01939985, lire en ligne).
  4. Yudell Luke, « On the computation of Log Z and Arc tan Z », Mathematical Tables and Other Aids to Computation, vol. 11,‎ , p. 16-18
  5. David K. Kahaner, « An Numerical quadrature by the ε-algorithm », Mathematics of computation, vol. 36,‎ , p. 689-693

Voir aussi[modifier | modifier le code]

Bibliographie[modifier | modifier le code]

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

Articles connexes[modifier | modifier le code]