# 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:
$~P(x)=a_{0}+a_{1}.x+a_{2}.x^{2}+\dotsb +a_{p}.x^{p}$
dont les racines λ$_k$ sont telles que |λ$_1$|> |λ$_2$|≥...≥|λ$_p$|
la méthode de Bernoulli génère la suite y$_n$ définie par les termes initiaux y$_0$, y$_1$...y$_{p-1}$ arbitraires, non tous nuls. Les termes suivants sont générés par récurrence:
$y_{n+p}={\frac {-1}{a_{p}}}\left(a_{p-1}.y_{n+p-1}+\dotsb +a_{1}.y_{n+1}+a_{0}.y_{n}\right){\text{ }}n=0,1\dotsb$
Le rapport des termes consécutifs y$_{n+1}$/y$_n$ 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
$y_{n+p}={\frac {-1}{a_{p}}}\left(a_{p-1}.y_{n+p-1}+\dotsb +a_{1}.y_{n+1}+a_{0}.y_{n}\right){\text{ }}n=0,1\dotsb$
est une suite récurrente linéaire à coefficients constants, dont la solution analytique est connue, et peut s'exprimer par:
$y_{n}=(c_{1,0}+c_{1,1}n+\dotsb +c_{1,k_{1}}n^{k_{1}})\lambda _{1}^{n}+(c_{2,0}+c_{2,1}n+\dotsb +c_{2,k_{2}}n^{k_{2}})\lambda _{2}^{n}+\dotsb +(c_{m,0}+c_{m,1}n+\dotsb +c_{m,k_{m}}n^{k_{m}})\lambda _{m}^{n}{\text{ }}n=0,1\dotsb$
où les λ$_i$ sont les m racines du polynôme
$~P(x)=a_{0}+a_{1}.x+a_{2}.x^{2}+\dotsb +a_{p}.x^{p}$
et k$_i$ la multiplicité de la racine λ$_i$ (avec k$_1$+k$_2$+..+k$_m$=p le degré du polynôme). Les c$_{i,j}$ sont des coefficients arbitraires, qui peuvent être déterminés grâce aux termes initiaux de la suite y$_0$, y$_1$...y$_{p-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 $|\lambda _{1}|>|\lambda _{2}|\geqslant |\lambda _{3}|\dotsb \geqslant |\lambda _{m}|$, la suite y$_n$ se comporte asymptotiquement comme:
$y_{n}=\lambda _{1}^{n}.c_{1,0}.[1+({\frac {\lambda _{2}}{\lambda _{1}}})^{n}.n^{k_{2}'}.(c_{2,k_{2}'}+O(1/n))+..]$
Le rapport de deux termes consécutifs se comporte asymptotiquement comme:
    ${\frac {y_{n+1}}{y_{n}}}=\lambda _{1}.{\frac {1+({\frac {\lambda _{2}}{\lambda _{1}}})^{n+1}.(n+1)^{k_{2}'}.(c_{2,k_{2}'}+O(1/n))+\dotsb }{1+({\frac {\lambda _{2}}{\lambda _{1}}})^{n}.n^{k_{2}'}.(c_{2,k_{2}'}+O(1/n))+\dotsb }}$
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 y$_0$, y$_1$...y$_{p-1}$, le coefficient c$_{1,0}$ serait égal à 0, le ratio des y$_n$ convergera vers λ$_2$.
Lorsque la racine dominante est de multiplicité k$_1$, le comportement asymptotique est
    $y_{n}=\lambda _{1}^{n}.[n^{k_{1}'}.(c_{1,k_{1}'}+O(1/n))+({\frac {\lambda _{2}}{\lambda _{1}}})^{n}.n^{k_{2}'}.(c_{2,k_{2}'}+O(1/n))+..]$
Les ki' sont au plus égaux à la multiplicité de la iéme racine (inférieur si certains c$_{i,ki}$ sont nuls).
Dans ce cas, le rapport de deux termes consécutifs se comporte asymptotiquement comme:
${\frac {y_{n+1}}{y_{n}}}=\lambda _{1}.{\frac {(n+1)^{k_{1}'}}{n^{k_{1}'}}}.{\frac {1+(c_{1,k_{1}'}+O(1/n))+({\frac {\lambda _{2}}{\lambda _{1}}})^{n+1}.(n+1)^{k_{2}'}.(c_{2,k_{2}'}+O(1/n))+\dotsb }{1+(c_{1,k_{1}'}+O(1/n))+({\frac {\lambda _{2}}{\lambda _{1}}})^{n}.n^{k_{2}'}.(c_{2,k_{2}'}+O(1/n))+\dotsb }}$
et la convergence vers λ$_1$ est en général logarithmique, sauf au cas où c$_{1,1}$, c$_{1,2}$, etc. c$_{1,k1}$=0 et c$_{1,0}$≠0, où elle est linéaire. Un certain choix judicieux de valeurs initiales y$_0$, y$_1$...y$_{p-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, y$_0$ = y$_1$...y$_{p-2}$ = 0 et y$_{p-1}$=1 garanti un coefficient c$_{1,0}$ non nul, donc une convergence vers la racine dominante.
Si une estimation de la racine dominante λ est disponible, l'initialisation y$_0$=1; y$_1$=λ; y$_2$=λ²...y$_{p-1}$=λ$^{p-1}$ permet d'améliorer la convergence en diminuant les valeurs des coefficients c$_{i,k}$ par rapport à c$_{1,0}$.
Le choix:
génère la suite y$_n$ 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:
$P_{n}=y_{n}.y_{n-2}-y_{n-1}^{2}$
$Q_{n}=y_{n}.y_{n-3}-y_{n-1}.y_{n-2}$
vérifient:
$\lim _{n\to \infty }{\frac {Q_{n}}{P_{n-1}}}=S$
$\lim _{n\to \infty }{\frac {P_{n}}{P_{n-1}}}=P$
où les deux racines dominantes complexes conjugués sont racines de:
$\lambda ^{2}-S\lambda +P=0$
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.
Dans le cas particulier ou les deux racines dominantes sont réelles et de signe opposé (λ$_2$ = -λ$_1$), et λ$_3$ la racine dominante suivante, on a:
$\lim _{n\to \infty }{\frac {y_{2n}}{y_{2n-2}}}=\lambda _{1}^{2}$ et $\lim _{n\to \infty }{\frac {y_{2n+1}}{y_{2n-1}}}=\lambda _{3}^{2}$

### 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 x$^p$.P(1/x), obtenu en reversant l'ordre des coefficients:
$~x^{p}.P(x^{-1})=a_{p}+a_{p-1}.x+a_{p-2}.x^{2}+\dotsb +a_{0}.x^{p}$
En appliquant la méthode de Bernoulli avec les coefficients du polynôme d'origine écrit à rebours,
    $y_{n+p}={\frac {-1}{a_{0}}}\left(a_{1}.y_{n+p-1}+\dotsb +a_{p-1}.y_{n+1}+a_{p}.y_{n}\right){\text{ }}n=0,1\dotsb$
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 y$_n$ 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 (Méthode de Durand-Kerner (en)).
Par exemple, la méthode d'Aitken consiste à déduire de la suite y$_k$, 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:
    $\lim _{k\to \infty }D_{m}^{(k)}=\lambda _{1}.\lambda _{2}.\lambda _{3}.\lambda _{4}....\lambda _{m}.$
d'où
    $\lim _{k\to \infty }{\frac {D_{m}^{(k)}}{D_{m-1}^{(k)}}}=\lambda _{m}$
    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 y$_{n+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 y$_{n+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 z$_m$=y$_{m-1}$/y$_m$, la relation de récurrence de Bernoulli devient:
    ${\frac {-a_{p}}{z_{n+p}}}=a_{p-1}+a_{p-2}.z_{n+p-1}+a_{p-3}.z_{n+p-1}.z_{n+p-2}+a_{p-4}.z_{n+p-1}.z_{n+p-2}.z_{n+p-3}+\dotsb +a_{1}.z_{n+p-1}.z_{n+p-2}\dotsb .z_{n+2}+a_{0}.z_{n+p-1}.z_{n+p-2}\dotsb .z_{n+1}{\text{ }}n=1,2\dotsb$
La suite des z$_m$ se calculent facilement à partir de la fin, en emboîtant les multiplications à la manière de Horner.
Les z$_m$ initiaux correspondant à y$_0$ = y$_1$...y$_{p-2}$ = 0 et y$_{p-1}$=1 sont z1 = z2...z$_{p-1}$ = 1 et Z$_p$=0
ceux correspondant aux identités de Newton Girard sont:
Les estimations de la racine dominante se déduisent de la définition de z$_m$ :
    $\lim _{n\to \infty }{\frac {1}{z_{n}}}={\lambda _{1}}$
    $\lim _{n\to \infty }{\frac {1}{z_{n}.z_{n-1}}}.{\frac {z_{n}-z_{n-1}}{z_{n-1}-z_{n-2}}}=P$
    $\lim _{n\to \infty }{\frac {1}{z_{n}}}.{\frac {z_{n}-z_{n-2}}{z_{n-1}-z_{n-2}}}=S$
    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 y$_n$/y$_{n-1}$ ou z$_n$. 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 y$_0$=1, y$_1$=λ,y$_1$=λ² ...y$_{p-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 c$_{i,k}$ induite par le changement de valeurs initiales, qui tendent vers les coefficients idéaux pour la convergence (c$_{i,k}$=0 sauf c$_{1,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, Bairstow (en) 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.
Soit à résoudre l'équation f(x)=0, et soit x$_0$ une valeur estimée proche de la racine cherchée. Le but est de trouver une suite de valeur x$_n$ convergeant vers la racine de f(x). Les dérivées successives de f(x) jusqu'à l'ordre p, sont calculées en x$_n$, et avec les notations: $a_{k}={\frac {f^{(k)}(x_{n})}{k!}}{\text{ }}k=0,1\dotsb p$
La fonction peut être approximée localement par son développement de Taylor, autour de x$_n$ :
$f(x_{n+1})=f(x_{n}+h)=a_{0}+a_{1}h+a_{2}h^{2}+\cdots +a_{p}h^{p}=0$
L'approximation suivante x$_{n+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 x$_{n+1}$ = x$_n$ +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é: R$_k$=a$_0$$^k$.y$_k$.
L'itération de la fonction se calcule par: $x_{n+1}=x_{n}+a_{0}.{\frac {R_{p-1}}{R_{p}}}$ La méthode est d'ordre p+1 pour les fonctions C$^{p+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:$x_{n+1}=x_{n}+a_{0}\cdot {\frac {R_{1}}{R_{2}}}=x_{n}+a_{0}\cdot {\frac {a_{1}}{-a_{1}^{2}+a_{2}a_{0}}}=x_{n}-{\frac {f_{n}f'_{n}}{{f'_{n}}^{2}-{\frac {f_{n}f''_{n}}{2}}}}$ , est attribuée à Halley. Elle est d'ordre 3 pour les racine simples, et d'ordre 1 pour les racines multiples.
La formule suivante, d'ordre 4, s'explicite en $x_{n+1}=x_{n}+a_{0}\cdot {\frac {R_{2}}{R_{3}}}=x_{n}-a_{0}\cdot {\frac {a_{1}^{2}-a_{2}a_{0}}{a_{1}^{3}-2a_{0}a_{1}a_{2}+a_{0}^{2}a_{3}}}=x_{n}-f_{n}\cdot {\frac {{f'_{n}}^{2}-{\frac {f_{n}f''_{n}}{2}}}{{f'_{n}}^{3}-f_{n}f'_{n}f''_{n}+{\frac {f_{n}^{2}f'''_{n}}{6}}}}$, correspond à la méthode de Householder d'ordre 4.
Il est possible d'exprimer directement les ratio R$_{p-1}$/R$_p$=w$_p$ à l'aide de:
${(w_{p})}^{-1}=a_{1}+a_{2}.a_{0}.w_{p-1}+a_{3}.a_{0}^{2}.w_{p-1}.w_{p-2}+a_{4}.a_{0}^{3}.w_{p-1}.w_{p-2}.w_{p-3}+\dotsb +a_{p}.a_{0}^{p-1}w_{p-1}.w_{p-2}\dotsb .w_{2}.w_{1}=\sum _{k=1}^{p}{\bigla_{k}\prod _{i=1}^{k-1}a_{0}w_{p-i}{\bigr )}{\text{ ; }}p=1,2\dotsb$
avec $x_{n+1}=x_{n}+a_{0}.w_{p}$
Les w$_k$ peuvent se calculer par la boucle suivante:
Avec l'initialisation par les identités de Newton Girard, en changeant la notation R$_k$ en T$_k$ :
L'itération de la fonction se calcule par: $x_{n+1}=x_{n}+a_{0}.{\frac {T_{p-1}}{T_{p}}}$ La méthode est d'ordre p pour les fonctions C$^{p+1}$, et reste d'ordre p pour des racines multiples.
Les T$_k$ peuvent se calculer par la boucle suivante:
En explicitant les formules, pour p=2, on obtient:
$x_{n+1}=x_{n}+a_{0}\cdot {\frac {T_{1}}{T_{2}}}=x_{n}+a_{0}\cdot {\frac {a_{1}}{2a_{2}a_{0}-a_{1}^{2}}}=x_{n}-{\frac {f_{n}f'_{n}}{(f'_{n})^{2}-f''_{n}f_{n}}}$, attribuée à Schröder, d'ordre 2, même pour les racines multiples.
La formule suivante, s'explicite en $x_{n+1}=x_{n}+a_{0}\cdot {\frac {T_{2}}{T_{3}}}=x_{n}-a_{0}\cdot {\frac {a_{1}^{2}-2a_{2}a_{0}}{a_{1}^{3}-3a_{0}a_{1}a_{2}+3a_{0}^{2}a_{3}}}=x_{n}-2f_{n}\cdot {\frac {{f'_{n}}^{2}-f_{n}f''_{n}}{2{f'_{n}}^{3}-3{f_{n}f'_{n}f''_{n}}+{f_{n}^{2}f'''_{n}}}}$ , est d'ordre 3, y compris pour les racines multiples.
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.

### Extension du domaine de convergence

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 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 3e 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 x$_{n+1}$ sera celle ayant la valeur de la fonction la plus proche de zéro.
- 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. 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 (étendant le domaine de convergence de la méthode par la possibilité de prendre en compte un pôle à proximité).
Les coefficients de l'approximant de Padé [p-1;1], connaissant le développement de Taylor à l'ordre p, s'expriment:
Pour prendre en compte un pôle éventuel à proximité de l'origine du développement de Taylor d'ordre p, il y aura donc à remplacer dans la méthode de Bernoulli les coefficients a$_k$ par les coefficients b$_k$. Ceux ci correspondent au numérateur de l'approximant de Padé, qui est de degré p-1. Les coefficients c$_0$ et c$_1$ sont ceux du dénominateur de l'approximant de Padé, qui ne servent ici que d'intermédiaire de calcul.
De même, si un couple de pôle réels ou complexes sont suspectés à proximité, les coefficients de l'approximant de Padé [p-2;2] s'expriment:
Dans ce cas, les coefficients b$_k$ (de degré p-2) sont à utiliser en lieu et place des a$_k$ dans la méthode de Bernoulli.
Les coefficients c$_0$ , c$_1$ et c$_2$ sont ceux du dénominateur de l'approximant de Padé [p-2;2]: ses racines représentent la localisation supposée de 2 pôles de la fonction à approximer. Si certains pôles de cette fonction sont connus de manière explicite, il est possible de remplacer le trinôme au dénominateur par celui ayant pour racine les pôles connus de la fonction. Dans ce cas, les coefficients b$_k$ seront à calculer jusqu'à l'ordre p.

## Exemples

* Exemple 1: racine réelles d'un polynôme
Estimation des racines du polynôme $x^{4}-22x^{3}+149x^{2}-308x+180$ , dont les racines sont: 1, 2, 9 et 10.
la récurrence pour la racine dominante est donc: $y_{n+4}=22y_{n+3}-149y_{n+2}+308y_{n+1}-180y_{n}$
et la récurrence pour la racine dominée est: $y_{n+4}={\frac {1}{180}}\left(308y_{n+3}-149y_{n+2}+22y_{n+1}-y_{n}\right)$
les valeurs initiales sont: y$_0$ = y$_1$ = y$_2$ = 0 et y$_3$ = 1
| | Racine dominée | Racine dominée | Racine dominée | Racine dominée | Racine dominée | Racine dominée | Racine dominante | Racine dominante | Racine dominante | Racine dominante | Racine dominante | Racine dominante |
| n | yn | Pn | Qn | λ$_4$ | S | P | yn | Pn | Qn | λ$_1$ | S | P |
| 3 | 1 | - | - | - | - | - | 1 | - | - | - | - | - |
| 4 | 1,71111111111111 | - | - | 0,58441558442 | - | - | 22 | - | - | 22 | - | - |
| 5 | 2,100123456790 | -1 | - | 0,81476691553 | - | - | 335 | -1 | - | 15,2272727273 | - | - |
| 6 | 2,299347050754 | -0,8277777778 | -1,2941975309 | 0,91335644878 | 1,5634601044 | 0,8277777778 | 4400 | -149 | -2970 | 13,1343283582 | 19,9328859 | 149 |
| 7 | 2,399583005639 | -0,4760802469 | -0,7229595336 | 0,95822776097 | 1,5185665406 | 0,5751304996 | 53481 | -15425 | -297418 | 12,1547727273 | 19,2815559 | 103,52348993 |
| 8 | 2,449780333960 | -0,2475763032 | -0,3726329637 | 0,97950945739 | 1,5051237092 | 0,5200306141 | 620202 | -1443865 | -27548730 | 11,5966791945 | 19,0798516 | 93,605510535 |
| 9 | 2,474888814884 | -0,1251034151 | -0,1878229595 | 0,98985470346 | 1,5013415850 | 0,5053125581 | 6970675 | -131328561 | -2498053162 | 11,2393623368 | 19,0214005 | 90,956260454 |
| 10 | 2,487444246098 | -0,0627225436 | -0,0941050070 | 0,99495247733 | 1,5003378627 | 0,5013655589 | 76624900 | -11851851129 | -225250299450 | 10,9924648617 | 19,0054952 | 90,245800599 |
| 11 | 2,493722104011 | -0,0313826501 | -0,0470765736 | 0,99748253508 | 1,5000827961 | 0,5003408395 | 828512861 | -1067393725825 | -20281941389578 | 10,8125799968 | 19,0013684 | 90,061351109 |
| 12 | 2,496861049779 | -0,0156939348 | -0,0235412146 | 0,99874284323 | 1,5000199039 | 0,5000831586 | 8845504382 | -96081412658825 | -1825578864841048 | 10,6763633956 | 19,0003333 | 90,014968548 |
| 13 | 2,498430524631 | -0,0078472805 | -0,0117709577 | 0,99937181569 | 1,5000047187 | 0,5000199474 | 93498427815 | -8647672122093568 | -1,643064610E+017 | 10,5701635291 | 19,0000799 | 90,003590526 |
| 14 | 2,499215262285 | -0,0039236773 | -0,0058855203 | 0,99968600638 | 1,5000011070 | 0,5000047238 | 980374738200 | -7,782978439E+017 | -1,478767375E+019 | 10,4854676288 | 19,0000189 | 90,000850285 |
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 :
Accélération de la convergence de λ1 par le Δ² :
| λ | Delta-2 1 | Delta-2 2 | Delta-2 3 |
| 22 | | | |
| 15,2272727273 | | | |
| 13,1343283582 | 12,19829858457 | | |
| 12,1547727273 | 11,29296300957 | | |
| 11,5966791945 | 10,85766044827 | 10,45452212745 | |
| 11,2393623368 | 10,60345511933 | 10,24662830085 | |
| 10,9924648617 | 10,44040269890 | 10,14873793204 | 10,06162681301 |
| 10,8125799968 | 10,32970722817 | 10,09566977349 | 10,03283865839 |
| 10,6763633956 | 10,25145613172 | 10,06272590565 | 10,00879613272 |
| 10,5701635291 | 10,19442606828 | 10,04116170382 | 10,00029804399 |
| 10,4854676288 | 10,15188286476 | 10,02694729028 | 9,999456763456 |
La nouvelle estimation de la plus grande racine peut être utilisé pour réinitialiser une nouvelle suite de Bernoulli avec comme valeur initiales y$_0$=1, y$_1$=λ,y$_1$=λ² ...y$_{p-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.
itération de l'accélération de la suite précédente :
| Suite n°2 | Suite n°3 |
| 9,99949585659666 | 10,0000000712436 |
| 9,99954277270112 | 10,0000000646106 |
| 9,99958773806757 | 10,000000058254 |
| 9,99962879690367 | 10,0000000524501 |
| 9,99966587396396 | 10,0000000472094 |
| 9,99969927030643 | 10,0000000424893 |
| 9,99972933388907 | 10,0000000382406 |
| 10,0000000767712 | 10,000000000001 |
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 7e 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:
$x^{4}-18,1670889561217x^{3}+91,2661704273649x^{2}-79,5299757817653x+3,16421413664943$
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:
| n | yn | λ$_4$ |
| 3 | 1 | - |
| 4 | 25,1341952052522 | 0,998014195022212 |
| 5 | 602,884534091775 | 0,999917659754444 |
| 6 | 14433,807501309 | 0,999996679810777 |
| 7 | 345536,985274708 | 0,999999866753945 |
| 8 | 8271929,81501342 | 0,99999999465651 |
| 9 | 198024574,44086 | 0,999999999785737 |
| 10 | 4740578409,91793 | 0,999999999991409 |
| 11 | 113486337337,666 | 0,999999999999656 |
| 12 | 2716788469371,51 | 0,999999999999986 |
| 13 | 65038133756546,3 | 1 |
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:
$a_{0}=x-S-2;a_{1}=1-C;a_{2}=S/2;a_{3}=C/6;...a_{n}=-a_{n-2}/(n.(n-1))$
Partant de x$_0$=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.
Zéro de la fonction f(x)=x-2sin(x)-x à l'aide de la méthode Rp d'ordre 8, en partant de x$_0$=0 :
| n | x$_n$ | a$_0$.R$_0$/R$_1$ | a$_0$.R$_1$/R$_2$ | a$_0$.R$_2$/R$_3$ | a$_0$.R$_3$/R$_4$ | a$_0$.R$_4$/R$_5$ | a$_0$.R$_5$/R$_6$ | a$_0$.R$_6$/R$_7$ | a$_0$.R$_7$/R$_8$ | a$_0$.R$_8$/R$_9$ | Description |
| 0 | 0 | -2 | -2 | 6,00000000000002 | -0,399999999999999 | -1,21951219512195 | -2,7032967032967 | 9,16546762589932 | -0,222222222222222 | -1,28154345227932 | non convergence vers λ$_1$ |
| | _ | | | 2,35294117647059 | 1,85294117647059 | 3,95647615090906 | 2,32778993119059 | 2,80500672889401 | 2,86233149074198 | 2,62795658279164 | recherche de λ$_3$ |
| 1 | 2,80500672889401 | -0,050029406074024 | -0,050317307537598 | -0,050332858595305 | -0,050332971748262 | -0,050332974623125 | -0,050332974647328 | -0,050332974647718 | | | convergence vers λ$_1$ |
| 2 | 2,75467375424629 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | | | convergence vers λ$_1$ |
| 3 | 2,75467375424629 | | | | | | | | | | convergence atteinte |
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 x$_1$ 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.