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