Méthode de Bernoulli

Cet article est une ébauche concernant l’analyse.

Vous pouvez partager vos connaissances en l’améliorant (comment ?) selon les recommandations des projets correspondants.

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 ^[1]. 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 modifier

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 |λ₁|> |λ₂|≥...≥|λ_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:

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

Propriétés modifier

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₁+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 $|\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 λ₁.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

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 λ₁ 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).

Variantes modifier

Choix des valeurs initiales modifier

Le choix le plus simple, y₀ = y₁...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₀=1; y₁=λ; y₂=λ²...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:

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

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 modifier

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.

Racine de plus petit module modifier

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 modifier

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-λ₁) à 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é^[2] 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 ^[3] , l'algorithme q-d de H.Rutishauser^[4] 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:

 $d_{0}^{(k)}=1$ 
 $d_{1}^{(k)}=y_{k}$ 
 $d_{m+1}^{(k)}={\frac {(d_{m}^{(k)})^{2}-d_{m}^{(k-1)}.d_{m}^{(k+1)}}{d_{m-1}^{(k)}}}$ 
 $D_{m}^{(k)}={\frac {d_{m}^{(k+1)}}{d_{m}^{(k)}}}$

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 modifier

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 modifier

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^[5] 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₀ = y₁...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:

 ${\frac {-a_{p}}{z_{1}}}={\frac {a_{p-1}}{p}}$ 
 ${\frac {-a_{p}}{z_{2}}}=a_{p-1}-2{\frac {a_{p-2}}{a_{p-1}}}$ 
 ${\frac {-a_{p}}{z_{3}}}=a_{p-1}+\left(a_{p-2}-3{\frac {a_{p-3}}{a_{p-1}}}\right).z_{2}$ 
 ${\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}$ 
 ${\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}$ 
… 
 ${\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}$

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 modifier

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₀=1, y₁=λ,y₁=λ² ...y_p-1=λ^p-1, où λ est l'estimation de λ₁ 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 modifier

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: $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₀^k.y_k.

 $R_{0}=1$ 
 $R_{1}=-a_{1}$ 
 $R_{2}=-\left(a_{1}.R_{1}+a_{2}.a_{0}.R_{0}\right)$ 
 $R_{3}=-\left(a_{1}.R_{2}+a_{2}.a_{0}.R_{1}+a_{3}.a_{0}^{2}.R_{0}\right)$ 
 $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)$ 
… 
 $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: $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.

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}{\bigl (}a_{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:

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 :

 $T_{1}=-a_{1}$ 
 $T_{2}=-\left(a_{1}.T_{1}+2.a_{2}.a_{0}\right)$ 
 $T_{3}=-\left(a_{1}.T_{2}+a_{2}.a_{0}.T_{1}+3.a_{3}.a_{0}^{2}\right)$ 
 $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)$ 
… 
 $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: $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:

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

Exemples modifier

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₀ = y₁ = y₂ = 0 et y₃ = 1


	Racine dominée						Racine dominante
n	yn	Pn	Qn	λ₄	S	P	yn	Pn	Qn	λ₁	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₀=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.

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

$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	λ₄
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, 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.

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
n	x_n	a₀.R₀/R₁	a₀.R₁/R₂	a₀.R₂/R₃	a₀.R₃/R₄	a₀.R₄/R₅	a₀.R₅/R₆	a₀.R₆/R₇	a₀.R₇/R₈	a₀.R₈/R₉	Description
0	0	-2	-2	6,00000000000002	-0,399999999999999	-1,21951219512195	-2,7032967032967	9,16546762589932	-0,222222222222222	-1,28154345227932	non convergence vers λ₁
	_			2,35294117647059	1,85294117647059	3,95647615090906	2,32778993119059	2,80500672889401	2,86233149074198	2,62795658279164	recherche de λ₃
1	2,80500672889401	-0,050029406074024	-0,050317307537598	-0,050332858595305	-0,050332971748262	-0,050332974623125	-0,050332974647328	-0,050332974647718			convergence vers λ₁
2	2,75467375424629	0	0	0	0	0	0	0			convergence vers λ₁
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 (λ₁ 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.

Notes modifier

↑ Daniel Bernoulli, « Observationes de seriebus recurrentibus », Commentarii Academiae Scientarium Imperialis Petropolitanae, Petropolis, vol. 3,‎ 1728, p. 85-100[1]
↑ J. H. Wilkinson, « The evaluation of zeros of ill-conditionned polynomials », Numerische Mathematik, vol. 1,‎ 1959, p. 150-166
↑ Alexander Craig Aitken, « On Bernoulli’s Numerical Solution of Algebraic Equations », Proceedings of the Royal Society of Edinburgh, vol. 46,‎ 1927, p. 289-305
↑ Henrici, Watkins, « Finding zeros of a polynomial by the Q-D algorithm », Communications of the ACM, vol. 8, issue 9,‎ 1965, p. 570-574
↑ Herbert S. Wilf, « The numercical solution of polynomial equations », Mathematical methods for digital computers,‎ 1960, p. 233-241
↑ E. Schröder, « Üeber unendlich viele Algorithmen zur Auflosung der Gleichungen », Mathematische Annalen,‎ 1870, p. 317-365

Références modifier

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 modifier

Portail de l'analyse

[1] Daniel Bernoulli, « Observationes de seriebus recurrentibus », Commentarii Academiae Scientarium Imperialis Petropolitanae, Petropolis, vol. 3,‎ 1728, p. 85-100[1]

[2] J. H. Wilkinson, « The evaluation of zeros of ill-conditionned polynomials », Numerische Mathematik, vol. 1,‎ 1959, p. 150-166

[3] Alexander Craig Aitken, « On Bernoulli’s Numerical Solution of Algebraic Equations », Proceedings of the Royal Society of Edinburgh, vol. 46,‎ 1927, p. 289-305

[4] Henrici, Watkins, « Finding zeros of a polynomial by the Q-D algorithm », Communications of the ACM, vol. 8, issue 9,‎ 1965, p. 570-574

[5] Herbert S. Wilf, « The numercical solution of polynomial equations », Mathematical methods for digital computers,‎ 1960, p. 233-241

[6] E. Schröder, « Üeber unendlich viele Algorithmen zur Auflosung der Gleichungen », Mathematische Annalen,‎ 1870, p. 317-365

[1]

[2]

[3]

[4]

[5]

[6]