Méthode de Bernoulli

Soit le polynôme:

${\displaystyle ~P(x)=a_{0}+a_{1}.x+a_{2}.x^{2}+\dotsb +a_{p}.x^{p}}$

dont les racines λ_k sont telles que |λ₁|> |λ₂|≥...≥|λ_p|

la méthode de Bernoulli génère la suite y_n définie par les termes initiaux y₀, y₁...y_p-1 arbitraires, non tous nuls. Les termes suivants sont générés par récurrence:

${\displaystyle 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 λ₁, 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.

la suite

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

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

${\displaystyle ~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₁+k₂+..+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₀, y₁...y_p-1.

Les racines ayant été numérotées suivant leur module décroissant. Si le polynôme possède une racine dominante λ₁ unique, c'est-à-dire tel que ${\displaystyle |\lambda _{1}|>|\lambda _{2}|\geqslant |\lambda _{3}|\dotsb \geqslant |\lambda _{m}|}$, la suite y_n se comporte asymptotiquement comme:

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

${\displaystyle {\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 λ₁.La convergence est linéaire et dépend du ratio entre les 2 racines dominantes (la convergence sera d'autant plus rapide que λ₁ sera grand par rapport à λ₂). Si, par un choix particulier de valeurs initiales y₀, y₁...y_p-1, le coefficient c_1,0 serait égal à 0, le ratio des y_n convergera vers λ₂.

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

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

${\displaystyle {\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 λ₁ 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₀, y₁...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).

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é.

${\displaystyle {\frac {-a_{p}}{z_{1}}}={\frac {a_{p-1}}{p}}}$
${\displaystyle {\frac {-a_{p}}{z_{2}}}=a_{p-1}-2{\frac {a_{p-2}}{a_{p-1}}}}$
${\displaystyle {\frac {-a_{p}}{z_{3}}}=a_{p-1}+\left(a_{p-2}-3{\frac {a_{p-3}}{a_{p-1}}}\right).z_{2}}$
${\displaystyle {\frac {-a_{p}}{z_{4}}}=a_{p-1}+a_{p-2}.z_{3}+\left(a_{p-3}-4{\frac {a_{p-4}}{a_{p-1}}}\right).z_{3}z_{2}}$
${\displaystyle {\frac {-a_{p}}{z_{5}}}=a_{p-1}+a_{p-2}.z_{4}+a_{p-3}.z_{4}z_{3}+\left(a_{p-4}-5{\frac {a_{p-5}}{a_{p-1}}}\right).z_{4}z_{3}z_{2}}$
… 
${\displaystyle {\frac {-a_{p}}{z_{p}}}=a_{p-1}+a_{p-2}.z_{p-1}+a_{p-3}.z_{p-1}z_{p-2}+\dotsb +\left(a_{1}-p.{\frac {a_{0}}{a_{p-1}}}\right).z_{p-1}z_{p-2}\dotsb z_{3}z_{2}}$

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 x₀ 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: ${\displaystyle 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 :

${\displaystyle 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₀^k.y_k.

${\displaystyle R_{0}=1}$
${\displaystyle R_{1}=-a_{1}}$
${\displaystyle R_{2}=-\left(a_{1}.R_{1}+a_{2}.a_{0}.R_{0}\right)}$
${\displaystyle R_{3}=-\left(a_{1}.R_{2}+a_{2}.a_{0}.R_{1}+a_{3}.a_{0}^{2}.R_{0}\right)}$
${\displaystyle R_{4}=-\left(a_{1}.R_{3}+a_{2}.a_{0}.R_{2}+a_{3}.a_{0}^{2}.R_{1}+a_{4}.a_{0}^{3}.R_{0}\right)}$
… 
${\displaystyle R_{p}=-\left(a_{1}.R_{p-1}+a_{2}.a_{0}.R_{p-2}+\dotsb +a_{p-1}.a_{0}^{p-2}.R_{1}+a_{p}.a_{0}^{p-1}.R_{0}\right)=-\sum _{k=1}^{p}a_{k}a_{0}^{k-1}R_{p-k}}$

L'itération de la fonction se calcule par: ${\displaystyle 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:${\displaystyle 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.

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

${\displaystyle {(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}{\bigl (}a_{k}\prod _{i=1}^{k-1}a_{0}w_{p-i}{\bigr )}{\text{$ ; }}p=1,2\dotsb }

avec ${\displaystyle x_{n+1}=x_{n}+a_{0}.w_{p}}$

Les w_k 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 R_k en T_k :

${\displaystyle T_{1}=-a_{1}}$
${\displaystyle T_{2}=-\left(a_{1}.T_{1}+2.a_{2}.a_{0}\right)}$
${\displaystyle T_{3}=-\left(a_{1}.T_{2}+a_{2}.a_{0}.T_{1}+3.a_{3}.a_{0}^{2}\right)}$
${\displaystyle T_{4}=-\left(a_{1}.T_{3}+a_{2}.a_{0}.T_{2}+a_{3}.a_{0}^{2}.T_{1}+4.a_{4}.a_{0}^{3}\right)}$
… 
${\displaystyle T_{p}=-\left(a_{1}.T_{p-1}+a_{2}.a_{0}.T_{p-2}+\dotsb +a_{p-1}.a_{0}^{p-2}.T_{1}+p.a_{p}.a_{0}^{p-1}\right)=-{\Bigl (}p.a_{p}.a_{0}^{p-1}+\sum _{k=1}^{p-1}a_{k}.a_{0}^{k-1}.T_{p-k}{\Bigr )}}$

L'itération de la fonction se calcule par: ${\displaystyle 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. Pour p=2, on obtient:

${\displaystyle 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[6], 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 3^e 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.

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.

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

Estimation des racines du polynôme ${\displaystyle 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: ${\displaystyle y_{n+4}=22y_{n+3}-149y_{n+2}+308y_{n+1}-180y_{n}}$

et la récurrence pour la racine dominée est: ${\displaystyle y_{n+4}={\frac {1}{180}}\left(308y_{n+3}-149y_{n+2}+22y_{n+1}-y_{n}\right)}$

les valeurs initiales sont: y₀ = y₁ = y₂ = 0 et y₃ = 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 y₀=1, y₁=λ,y₁=λ² ...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.

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 7^e 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:

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

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:

${\displaystyle 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, 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 λ₃, à 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 (λ₁ et λ₂), l'estimation de λ₃ 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₁ 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.

• 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.