Méthode de Bernoulli
En analyse numérique, la méthode de Bernoulli est une méthode de calcul numérique de la racine de plus grand module d'un polynôme, exposée par Daniel Bernoulli . C'est une méthode itérative ne nécessitant aucune estimation préalable de la racine, qui sous certaines conditions, est globalement convergente.

Présentation 
Soit le polynôme:

dont les racines λk sont telles que |λ1|> |λ2|≥...≥|λp| 

la méthode de Bernoulli génère la suite yn définie par les termes initiaux y0, y1...yp-1 arbitraires, non tous nuls. Les termes suivants sont générés par récurrence:

Le rapport des termes consécutifs yn+1/yn converge, sous certaines conditions, vers λ1, la racine "dominante" du polynôme P(x). 

Cette méthode peut être rattachée à celle de la puissance itérée pour estimer les valeurs propres des matrices, en l'appliquant à la matrice compagnon du polynôme.

Propriétés 
la suite   

est une suite récurrente linéaire à coefficients constants, dont la solution analytique est connue, et peut s'exprimer par:

où les λi sont les m racines du polynôme 

et ki la multiplicité de la racine λi (avec k1+k2+..+km=p le degré du polynôme). Les ci,j sont des coefficients arbitraires, qui peuvent être déterminés grâce aux termes initiaux de la suite y0, y1...yp-1.

Les racines ayant été numérotées suivant leur module décroissant. Si le polynôme possède une racine dominante λ1 unique, c'est-à-dire tel que
, 
la suite yn se comporte asymptotiquement comme:

Le rapport de deux termes consécutifs se comporte asymptotiquement comme:

et converge vers la racine dominante simple λ1.La convergence est linéaire et dépend du ratio entre les 2 racines dominantes (la convergence sera d'autant plus rapide que λ1 sera grand par rapport à λ2). Si, par un choix particulier de valeurs initiales y0, y1...yp-1, le coefficient c1,0 serait égal à 0, le ratio des yn convergera vers λ2.

Lorsque la racine dominante est de multiplicité k1, le comportement asymptotique est

Les ki' sont au plus égaux à la multiplicité de la iéme racine (inférieur si certains ci,ki sont nuls).  
Dans ce cas, le rapport de deux termes consécutifs se comporte asymptotiquement comme:

et la convergence vers λ1 est en général logarithmique, sauf au cas où c1,1, c1,2, etc. c1,k1=0 et c1,0≠0, où elle est linéaire. Un certain choix judicieux de valeurs initiales y0, y1...yp-1 permet d'obtenir ce comportement.

La convergence étant généralement assurée quel que soit les termes initiaux choisit qualifie la méthode de globalement convergente, et rend cette méthode quasiment insensible aux propagations d'erreur (équivalant à un léger changement des valeurs initiales).

Variantes

Choix des valeurs initiales 
Le choix le plus simple, y0 = y1...yp-2 = 0 et yp-1=1 garanti un coefficient c1,0 non nul, donc une convergence vers la racine dominante.

Si une estimation de la racine dominante λ est disponible, l'initialisation y0=1; y1=λ; y2=λ²...yp-1=λp-1 permet d'améliorer la convergence en diminuant les valeurs des coefficients ci,k par rapport à c1,0.

Le choix:
 
 
 
 
 … 
 

génère la suite yn de la somme des racines élevé à la puissance n+1, d'après les identités de Newton Girard. Avec ces valeurs initiales, la suite du rapport de deux termes consécutifs converge linéairement vers la racine dominante, même dans le cas d'une racine multiple.

Racines complexes ou de même module 

Pour les polynômes à coefficients réels, si les termes initiaux de la suite sont aussi réels, la suite générée ne peut quitter l'axe des réels. Au cas où la racine dominante est un couple de racines complexes conjugués, la suite des ratios successifs oscillent sans converger. Dans ce cas, les suites auxiliaires:

vérifient:

où les deux racines dominantes complexes conjugués sont racines de:

Cette variante s'applique aussi dans le cas de racines dominantes réelles. Lorsque les deux racines dominantes forment un couple isolé des autres, la convergence peut être plus rapide vers S et P que le vers la racine dominante seule. Les mêmes remarques s'appliquent lorsque la racine dominante est double.

Racine de plus petit module 

La racine de plus petit module de P(x) est l'inverse de la racine de plus grand module de xp.P(1/x), obtenu en reversant l'ordre des coefficients:

En appliquant la méthode de Bernoulli avec les coefficients du polynôme d'origine écrit à rebours, 

le ratio des termes consécutifs de la suite converge vers 1/λp, inverse de la racine de plus petit module. Ainsi, toutes les propriétés et méthodes spécifiques aux racines dominantes peuvent se transcrire aux racines de plus faible module. En cas de racine multiple ou complexe conjugué, les mêmes procédés que pour les racines dominantes s'appliquent.

Racines intermédiaires 
 Déflation du polynôme
Une fois déterminée avec suffisamment de précision la racine de plus grand module, la division du polynôme de départ par (x-λ1) à l'aide du schéma de Horner, permet d'obtenir un polynôme de degré n-1 avec les racines restantes à déterminer. La méthode de Bernoulli peut être appliquée à ce polynôme réduit, pour de proche en proche, déterminer toutes les racines. De la même manière, il est possible de retirer la racine de plus petit module, par la méthode de Bernoulli avec les coefficients écrits à rebours. En cas de racines dominantes complexes conjuguées, une déclinaison du schéma de Horner permet de retirer les deux racines simultanément en divisant le polynôme par x²-Sx+P.
Avec ce type de procédé, les polynômes successifs peuvent avoir des racines légèrement différentes du polynôme initial, à la suite de propagations d'erreur d'arrondi. Afin de limiter au mieux cette propagation d'erreur, il est conseillé  de supprimer les racines par ordre croissant des modules. Il est aussi conseillé de garder en mémoire le polynôme d'origine afin d'affiner les racines par un autre procédé. 

 Détermination de toutes les racines simultanément

Certains algorithmes complémentaires permettent, en se basant sur la suite yn générée par la méthode de Bernoulli, d'estimer toutes les racines du polynôme simultanément. Parmi ces méthodes, figure celle d'A. Aitken
 , l'algorithme q-d de H.Rutishauser ou l'ε-algorithme de P. Wynn. Contrairement à la méthode de Bernoulli, ces méthodes complémentaires sont très sensibles à la propagation des erreurs d'arrondi, et sauf précautions particulières, la précision obtenue pour certaines racines peut être médiocre. Pour limiter les instabilités numériques de ces méthodes, il est souhaitable de déterminer les plus grandes racines en partant de la variante de Bernoulli convergeant vers la racine dominante, et de déterminer les plus petites racines à l'aide de la variante convergeant vers la racine de plus petit module.
Les valeurs trouvées pourront servir de premières estimations pour des algorithmes de résolution nécessitant des estimations de toutes les racines simultanément ().

Par exemple, la méthode d'Aitken consiste à déduire de la suite yk, générée par la méthode de Bernoulli, des suites associées, dont les ratios convergent vers le produit des m plus grandes racines du polynôme:

 
 
 
 

vérifient:

d'où

Des variantes permettent de prendre en charge les racines complexes.

Mise en œuvre calculatoire 
Les principaux points de vigilance de la méthode de Bernoulli sont les risques de dépassement de capacité des nombres en virgule flottant lors de la génération de la suite, et la lenteur de convergence du procédé.

Gestion des limitations des nombres flottants 
La suite générée par la méthode de Bernoulli tend à suivre une progression géométrique dont la raison est la plus grande (ou la plus petite) des racines du polynôme. Si la valeur de cette racine est élevée, la suite générée peut rapidement dépasser la capacité des nombres en virgule flottante utilisés (overflow ou underflow). Afin d'éviter cet écueil, il existe plusieurs stratégies.

 normalisation des racines du polynôme
Il existe plusieurs formules permettant d'estimer la dimension de la plus grande (ou plus petite) racine du polynôme, à l'aide des coefficients (entre autres, les premiers termes de la méthode de Bernoulli).
Une fois cette estimation, λmax, connue, un changement de variable z=x/|λmax| permet de ramener la plus grande racine à une valeur voisine de l'unité. Pour éviter la propagation des erreurs de calcul lors de ce changement de variable, la valeur de λmax pourra être remplacée par la puissance de 2 la plus proche. De même, pour éviter que les produits intermédiaires lors du calcul de la suite de Bernoulli génèrent des dépassements de capacité, une normalisation de tous les coefficients par une puissance de 2 adéquate est à envisager, une fois le changement de variable effectué. Les coefficients d'origine seront à conserver comme référence.
 
 rectification périodique de la suite générée

Lors du calcul de la suite des yn+p, une vérification peut être faite en cas de risque prochain de dépassement de capacité (le ratio de deux termes consécutifs donnant la raison de la suite géométrique, une estimation du terme suivant est aisée). En cas de dépassement de capacité proche, une division des p derniers termes de la suite par une puissance de 2 adéquate permet de repartir avec des termes loin du dépassement de capacité. La convergence de la suite n'est pas altérée par ce procédé, car seuls les ratios interviennent dans l'estimation des racines.

 utilisation de la suite des ratios consécutifs
Seuls les ratios de termes  yn+p sont d'intérêt pratique, et non les termes individuellement. Une reformulation permet de calculer directement le ratio à l'aide d'une relation de récurrence. Ces ratios convergent directement vers la valeur des racines, évitant les risques de dépassement de capacité.

En introduisant zm=ym-1/ym, la relation de récurrence de Bernoulli devient:

La suite des zm se calculent facilement à partir de la fin, en emboîtant les multiplications à la manière de Horner. 

Les zm initiaux correspondant à y0 = y1...yp-2 = 0 et yp-1=1 sont z1 = z2...zp-1 = 1 et Zp=0

ceux correspondant aux identités de Newton Girard sont:

 
 
 
 
 
 … 
  

Les estimations de la racine dominante se déduisent de la définition de zm :

Où S et P sont la somme et le produit des deux racines dominantes. La même méthode, en renversant l'ordre des coefficients, peut être utilisée pour calculer l'inverse de la ou des racine(s) de plus petit module.

Accélération de la convergence 
La méthode de Bernoulli converge typiquement de manière linéaire vers la racine cherchée, il est utile d'accélérer ce procédé. Plusieurs méthodes peuvent être mises en œuvre:

 utilisation d'algorithme d'accélération  
Plusieurs algorithmes d'accélération de la convergence sont aptes à accélérer la convergence de la méthode de Bernoulli. Parmi ceux-ci certains sont particulièrement adaptés au vu de la forme du terme d'erreur de la méthode: le Delta-2 d'Aitken, et l'ε-algorithme. Le but de ces algorithmes est d'accélérer la convergence de la suite yn/yn-1 ou zn. Chaque utilisation d'un des deux algorithmes supprime le plus grand des termes du ratio des deux racines dominantes non encore supprimés. La répétition de l'utilisation de ces algorithmes permet d'obtenir une suite à convergence linéaire, mais avec un terme d'erreur de plus en plus petit. Les termes consécutifs de la suite de Bernoulli utilisés pour accélérer la convergence ne doivent pas subir de transformation (décalage du zéro, rectification périodique) sur la plage utilisée pour l'accélération.  Il est en revanche possible ensuite de réinitialiser la suite de Bernoulli à l'aide du résultat obtenu par l'accélération de la convergence. En prenant comme valeur initiale y0=1, y1=λ,y1=λ² ...yp-1=λp-1, où λ est l'estimation de  λ1 trouvé par l'accélérateur, la nouvelle suite générée converge plus rapidement que l'originale. Ceci est dû à la modification des coefficients ci,k induite par le changement de valeurs initiales, qui tendent vers les coefficients idéaux pour la convergence (ci,k=0 sauf c1,0).
 décalage de l'origine 
Le coefficient de linéarité de la convergence est proportionnel au rapport de la racine cherchée avec celle la plus proche. Si les deux racines dominantes sont voisines, la convergence sera d'autant plus lente. Lorsque la méthode est employée pour chercher la racine de plus faible module, une technique permet d'augmenter le rythme de convergence: le décalage de l'origine. En supposant que l'algorithme fournit à un certain stade une approximation grossière de la racine de plus petit module, le fait de décaler l'origine à cette valeur permet de diminuer le ratio entre les deux plus petites racines, et donc d'augmenter le rythme de convergence de l'algorithme par la suite. Ce décalage peut être renouvelé par la suite pour accélérer davantage la convergence. Le décalage systématique de l'origine toutes les p itérations permet d'obtenir une convergence d'ordre p. Le cumul des décalages successifs doit être additionné pour déterminer la racine du polynôme d'origine.  Les coefficients du polynôme avec le décalage d'origine est obtenu à l'aide du schéma d'Horner emboîté.
 bascule avec une autre méthode 
Le principal intérêt de la méthode de Bernoulli est de fournir une estimation d'une des racines du polynôme, sans estimation préalable. L’inconvénient est sa convergence seulement linéaire. D'autres méthodes, à l'inverse, ont une convergence rapide (typiquement d'ordre deux), mais nécessitent une bonne estimation de départ, sous peine de diverger ou entrer dans des cycles itératifs. La combinaison de la méthode de Bernoulli, pour fournir une valeur approchée d'une racine, avec une méthode à convergence rapide, pour polir cette estimation, est complémentaire. Parmi les méthodes à associer, figure celle de Newton ou Birge-Vieta,  ou Laguerre.

Extension à des fonctions non polynomiales 
La méthode de Bernoulli est basée sur les propriétés des polynômes, et leur lien avec les suites récurrentes linéaires. Mais il est possible d'étendre la méthode à des fonctions autres que polynomiales, en les approximant localement par un polynôme. Dans ce cas, la méthode recherchant la racine de plus petit module s'avère intéressante. La technique du décalage systématique de l'origine pourra être mise en œuvre pour obtenir des méthodes d'ordre quelconque; celui-ci prenant dans ce cas la forme d'une méthode de point fixe classique.

Si les dérivées successives de la fonction f(x) sont connues jusqu'à un certain ordre, il sera possible de remplacer la fonction par son développement de Taylor autour du point où la fonction est susceptible d'avoir un zéro. Il est aussi possible de remplacer le développement de Taylor par un des approximant de Padé correspondant, en cherchant la racine du numérateur par la méthode de Bernoulli.

Soit à résoudre l'équation f(x)=0, et soit x0 une valeur estimée proche de la racine cherchée. Le but est de trouver une suite de valeur xn convergeant vers la racine de f(x).
Les dérivées successives de f(x) jusqu'à l'ordre p, sont calculées en xn, et avec les notations:

La fonction peut être approximée localement par son développement de Taylor, autour de xn :

L'approximation suivante xn+1 est déterminée en appliquant la méthode de Bernoulli sur l'équation polynomiale en h, en recherchant la racine de plus petit module, puis xn+1 = xn +h.

Avec l'utilisation systématique du décalage d'abscisse, la méthode de Bernoulli classique risque de rencontrer des difficultés de calcul, à cause du coefficient constant du polynôme devenant très petit. Pour éviter cela, un changement de variable est opéré: Rk=a0k.yk.

 
 
 
 
 
 … 
 

L'itération de la fonction se calcule par:

La méthode est d'ordre p+1 pour les fonctions Cp+1, mais chute à l'ordre 1 pour des racines multiples. La première formule obtenue est celle de Newton Raphson. La formule suivante, s'écrit en l'explicitant: , est attribuée à Halley. Elle est d'ordre 3 pour les racine simples, et d'ordre 1 pour les racines multiples.

Il est possible d'exprimer directement les ratio Rp-1/Rp=wp à l'aide de:

avec 

Les wk peuvent se calculer par la boucle suivante:

FOR k=1 TO p
   temp=a(k)
   FOR i=1 TO k-1
      temp=temp*a(0)*w(i)+a(k-i)
   NEXT i
   w(k)=-1/temp
NEXT k
Avec l'initialisation par les identités de Newton Girard, en changeant la notation Rk en Tk :
 
 
 
 
 … 
 
L'itération de la fonction se calcule par:

La méthode est d'ordre p pour les fonctions Cp+1, et reste d'ordre p pour des racines multiples. Pour p=2, on obtient: 
 
, attribuée à Schröder, d'ordre 2, même pour les racines multiples.

Le fait de prolonger la suite de Bernoulli au delà de p, dérivée la plus élevée disponible, n'augmente pas l'ordre de convergence de la méthode, mais modifie à la marge le domaine de convergence de la méthode. La suite de Bernoulli convergeant vers la racine du polynôme d'approximation, et non la racine de la fonction elle-même. Cette prolongation peut être envisagée lorsque les ratios Rp/Rp+1 convergent lentement ou oscillent, a des fins de diagnostic sur les cas de divergence, ou d'extension du domaine de convergence. De même, utiliser les méthodes d'accélération, type Delta-2 dans ces cas de figure, peut s’avérer utile pour stabiliser les convergences oscillantes. 

La convergence globale dont bénéficie la méthode de Bernoulli sur les polynômes ne se généralise pas avec la version non polynomiale. Les causes pouvant générer des limitations de domaine de convergence sont: 

- présence d'un pôle dans un rayon plus proche que celui de la racine cherchée, par rapport au point de départ de l'itération. Le pôle limite le rayon de convergence de la série de Taylor, et la racine se trouve hors de cette limite. La convergence globale n'est plus garantie. Ce cas peut être détecté par des estimations erratiques des différents ratios Rp/Rp+1. Un palliatif est de retenir seulement ceux liés à un polynôme de degré impair, ou de remplacer le développement de Taylor par le numérateur d'un approximant de Padé équivalent (étendant le domaine de convergence par la possibilité de prendre en compte un pôle à proximité). 

- présence de racines complexes conjugués plus proche du point de départ des itérations que la racine réelle cherchée. Ce cas est détecté par l'oscillation des ratios Rp/Rp+1, et l'éventuelle convergence des ratios de la somme et du produit des plus faibles racines. Un contournement est de retenir comme première proposition R0/R1 (méthode de Newton), comme deuxième proposition d'évaluer la plus grande racine du développement de Taylor d'ordre 3 (dont les deux plus petites racines seraient complexes conjugués, et la plus grande, la seule racine réelle). La troisième proposition serait de calculer la  racine de plus petit module à l'aide des Rp déjà calculés, et de la méthode d'Aitken ou de l'epsilon algorithme. La proposition retenue pou xn+1 sera celle ayant la valeur de la fonction la plus proche de zéro.

L'usage de ce genre de méthode, d'ordre élevé, se justifie lorsque les dérivées d'ordre élevé sont de calcul simple, ou nécessite peu de calcul supplémentaire une fois les premières dérivées calculées. Ceci se présente aussi lorsque la fonction vérifie une équation différentielle simple.

Exemples 

 Exemple 1: racine réelles d'un polynôme

Estimation des racines du polynôme  , dont les racines sont: 1, 2, 9 et 10.

la récurrence pour la racine dominante est donc: 

et la récurrence pour la racine dominée est: 

les valeurs initiales sont: y0 = y1 = y2 = 0 et y3 = 1

La convergence vers la racine dominée "1" est linéaire mais relativement rapide: un chiffre exact est gagné toutes les deux ou trois itérations. Ce rythme s'explique par la racine la plus proche qui est le double de la racine dominée: l'erreur évolue comme une suite géométrique de rapport 1/2. La convergence vers S et P (la somme et le produit des deux plus faibles racines) est meilleure, du fait que les racines suivantes sont significativement éloignées. S et P convergent vers 1.5 et 0.5, dont l'inverse des racines de x²-Sx+P sont 1 et 2.

La convergence vers la racine dominante "10" est nettement plus lente que la racine dominée: le ratio avec la racine la plus proche valant 0.9, l'erreur est vaguement est multipliée par cette valeur à chaque itération. La convergence de S et P vers la somme et le produit des deux plus grandes racines est nettement plus rapide, les deux racines suivantes étant significativement éloignés. S et P convergent vers 19 et 90, dont les racines de x²-Sx+P sont 9 et 10.

La convergence particulièrement lente de la plus grande racine peut être accélérée par le Δ² d'Aitken, par exemple. Voici le résultat avec les données extraites du tableau précédent :

La nouvelle estimation de la plus grande racine peut être utilisé pour réinitialiser une nouvelle suite de Bernoulli  avec comme valeur initiales y0=1, y1=λ,y1=λ² ...yp-1=λp-1, qui convergera plus vite que la suit d'origine. Cette nouvelle suite pouvant à nouveau être accélérée.

Le tableau suivant continue le procédé précédent, en générant 7 termes de la suite de Bernoulli en initialisant les valeurs initiales comme indiqué. Le terme en gras est issu de l'application en cascade du Δ² de cette suite, qui servira de germe pour la prochaine suite. L'accélération de la convergence par rapport à la suite originale est nettement améliorée.

La convergence vers la plus petite racine peut être accélérée par un décalage de l'origine. Par exemple, en s'arrêtant à la  itération, une valeur approchée grossière de la racine est 0,95822776097. En décalant l'origine à cette abscisse, le polynôme décalé devient:

en recommençant les itérations avec ce polynôme, et en appliquant le décalage aux estimations trouvées par la méthode de Bernoulli, on trouve:

La convergence a été nettement accélérée. La nouvelle racine dominée étant de l'ordre de 0.05, et la racine la plus proche de l'ordre de 1, leur ratio est nettement plus petit que le ratio avant décalage, d'où l'accélération de la convergence. Ce procédé peut être répété en cumulant les décalages, déduites des nouvelles estimations de la racine dominée. Il est possible aussi de coupler le procédé au Δ² d'Aitken, en proposant comme décalage non pas l'estimation de Bernoulli mais celle du Δ².

 Exemple 2: zéro d'une fonction non polynomiale

La fonction f(x)=x-2*sin(x)-2 possède une racine réelle voisine de 2.75. Ses dérivées successives se calculent à très peu de frais, une fois évaluée S=2*sin(x) et C=2*cos(x), par:

Partant de x0=0, la présence de deux racines complexes conjugués à une distance d'environ 1.6 (donc plus proche de la racine réelle), devrait faire osciller les suites de Bernoulli. La détection de ces oscillations orientera le calcul vers l'estimation de λ3, à l'aide de la suite déjà calculée et de la méthode d'Aitken.

Les dérivées successives ont été calculées jusqu'à l'ordre 7, représentant une méthode théoriquement d'ordre 8. Dans le cas d'une suspicion de non convergence, la suite de Bernoulli a été poussée jusqu'à 9, sans ajout des dérivées 8 et 9. À la suite de la détection de non-convergence, et de suspicion de présence de deux racines complexes conjuguées proches (λ1 et λ2), l'estimation de λ3 utilise la suite déjà calculée et la méthode d'Aitken. L'itération suivante converge sans encombre, la racine réelle étant cette fois nettement plus proche de x1 que les deux racines complexes. À noter que sur cet exemple, les méthodes de Newton et Halley génèrent des itérations chaotiques, avant de converger après respectivement 17 et 77 itérations.

Notes

Références 
 Emile Durand, Solutions numérique des équations algébriques, Volume 1, Equations du type F(x)=0 , racines d'un polynôme,  éditions Masson, 1960.
 Francis B Hildebrand, Introduction to numerical analysis, second edition, McGraw-Hill, Inc, 1956.

Voir aussi 

Bernoulli
Bernoulli