Calcul numérique des racines d'un polynôme

I) Présentation

1) Objectif

On veut calculer (numériquement) toutes les racines complexes d'un polynôme à coefficients complexes.

2) Un problème difficile

Même un logiciel très performant comme Mathématica n'est pas au point lorsqu'il y a des racines multiples. Par exemple:

f [ x_ ] = Expand [ ( x - 1. ) 5 ( x - 2. ) 6 ] ; NSolve [ f [ x ] 0 ]

{{x -> 0.994927}, {x -> 0.997181}, {x -> 1.00011}, {x -> 1.00331}, {x -> 1.00385}, {x -> 1.97488}, {x -> 1.98512}, {x -> 1.99182}, {x -> 2.00666}, {x -> 2.01183}, {x -> 2.02496}}

Les solutions devraient être 1.0 (cinq fois) et 2.0 (6 fois)

f [ x_ ] = Expand [ ( x - 3.14 ) 5 ( x - 2.718 ) 6 ] ; NSolve [ f [ x ] 0 ]

{{x -> 2.90982}, {x -> 2.90982}, {x -> 2.90982}, {x -> 2.90982}, {x -> 2.90982}, {x -> 2.90982}, {x -> 2.90982}, {x -> 2.90982}, {x -> 2.90982}, {x -> 2.90982}, {x -> 2.90982}}

  Les solutions devraient être 3.14 (cinq fois) et 2.718 (6 fois)

  Aucun algorithme ne peut calculer correctement les n racines complexes d'un polynôme de degré n quelconque, en calculant uniquement sur les flottants usuels. Les problèmes liés aux calculs sur les nombres flottants (notamment les erreurs d'élminations) provoqueront toujours, pour des polynômes bien choisis, des résultats inaceptables.

  Nous allons étudier des idées qui donnent de bons résultats, mais qui auront aussi leurs limites, comme toutes les méthodes numériques.

3) Méthode

Le calcul des racines d'un polynôme p(x) de C[x] sera fait en trois phases:
     (1) Calcul du polynôme q(x) unitaire ayant les mêmes racines que le polynôme p(x), mais toutes d'ordre 1
     (2) Calcul des racines du polynôme q(x)
     (3) Calcul des racines du polynôme p(x) avec leur ordre de multiplicité

4) Mathématiques utiles

  Chacune des trois phases du calcul des racines d'un polynôme p(x) de C[x] nécessite l'utilisation d'un théorème de mathématique. On admet ces trois théorèmes, seul le second est délicat à prouver.

a) Phase (1) : calcul du polynôme q(x) ayant les mêmes racines que le polynôme p(x), mais toutes d'ordre 1

Théorème 1

Soit p C [ x ]   un polynôme non constant. Alors:
    le polynôme q ( x ) = p ( x ) pgcd ( p ( x ) , p ' ( x ) ) a les mêmes racines que le polynôme p(x), mais toutes d'ordre 1.

b) Phase 2: calcul des racines du polynôme q(x)

Théorème 2

Soit p C [ x ] un polynôme non constant. On définit la fonction m : C R + z p ( z ) . Alors:
    les complexes a m admet un minimum local sont les racines du polynôme p ( x )

Remarque 1: En une racine a de p, m(a) = 0. Comme m est une fonction à valeurs positive, les racines de p sont donc des minimums locaux de la fonction m. Mais il n'est pas évident c'est que ce sont les seuls minimums locaux.

Remarque 2: La version réelle du théorème est fausse: a = 1 est un minimum local (et même global) de la fonction m ( x ) = | 1 + ( x - 1 ) 2 | définie sur R, mais a = 1 n'est pas une racine de p ( x ) = 1 + ( x - 1 ) 2 .

c) Phase 3 : calcul des racines du polynôme p(x) avec leur ordre de multiplicité

Théorème 3

Soient p C [ x ] un polynôme non constant, a C et n N . Alors:
     a est une racine d'ordre n de p p ( a ) = p ' ( a ) = ... = p ( n - 1 ) ( a ) = 0 et p ( n ) ( a ) 0

5) Modélisation d'un polynôme

Le polynôme non nul de degré n N     p ( x ) = Σ k = 0 n a k x k de C[x] est la liste de complexes   p = [ a 0 , a 1 , ... , a n ] (avec a n 0 ) .
Le polynôme nul est la liste [0], et son degré sera -1 .

On suppose disposer des fonctions de base suivantes pour le calcul sur les polynômes

6) Fonctions de base

from polynomes import *

On dispose alors de toutes les fonctions de base écrites dans le fichier polynomes.py : degré, opérations usuelles, division euclidienne, dérivées, pgcd.

def deg(p):
    """ Degré du polynome p. De plus, deg(0) = -1"""

def valeur(p, x):
    """ p(x) calculé par l'algorithme de Horner"""

def somme(p, q):
    """ polynôme p+q"""

def fois(a, p):
    """ Polynôme a.p """

def difference(p,q):
    """ polynôme p-q"""

def produit(p, q):
    """ polynôme produit r = p x q """

def puissance(p,n):
    """Polynôme p^n"""

def composition(p, q):
    """ polynôme poq """

def derivee_1(p):
    """ Dérivée du polynôme p """

def derivee(p, n):
    """ Dérivée nième du polynôme p """

def norme(p):
    """ Plus grand module des coefficients de p"""

def nettoyer(p, m):
    """ Mise à zéro des coefficients de p tels que |c / m| < eps = 2^(-30)"""

def reste(a, b):
    """ reste de la division de a par b"""

def quotient(a, b):
    """ quotient de la division de a par b"""

def unitaire(p):
    """ Polynôme unitaire colinéaire à p """

def pgcd(p, q):
    """ PGCD unitaire des polynômes p et q"""

7) Nombres complexes

  Le complexe z = a + i b se note en Python a + b j , y compris si b = ± 1 . Les parties réelles et imaginaires de z sont z.real et z.imag, le module de z est abs(z).

II) Phase (1) : calcul du polynôme q(x) ayant les mêmes racines que le polynôme p(x), mais toutes d'ordre 1

Q1) En expolitant le théorème 1, écrire une fonction polynome_racines_simples(p) (avec p polynôme) calculant le polynôme q unitaire ayant les mêmes racines que p, mais toutes d'ordre 1.

  Par exemple, avec p = produit ( puissance ( [ - 4 , 3 + 5 j , 1 - 7 j ] , 5 ) , puissance ( [ 1 , - 2 + 1 j ] , 3 ) ) , vérifier que
q = polynome_racines _simples ( p ) = unitaire ( produit ( [ - 4 , 3 + 5 j , 1 - 7 j ] , [ 1 , - 2 + 1 j ] ) )

Q2) Les erreurs d'arrondis dans les calculs peuvent fausser le calcul du polynôme q(x).

  Ecrire une fonction pol1(n) (avec n entier) calculant le polynôme p1(n,x)=x(x-1)...(x-n) et vérifier que le polynôme q(x) n'est pas correctement calculé pour p = p1 ( 4 , x ) 3

  Modifier la constante eps=2**(-30) de la fonction nettoyer(p, m) et essayer avec eps = 2**(-26) : le calcul de q(x) est cette fois correct.

  L'ajustement de cette constante eps est délicat: que se passera-t-il si eps est trop grande (par exemple eps = 2**(-10) ) ?

III) Phase (2) : calcul des racines du polynôme q(x)

1) Calcul de toutes les racines

Le calcul des racines d'un polynôme p(x) scindé est résumé par l'algorithme suivant:

tant que deg(p) > 0:
    calculer une racine a de p(x)
    p(x) = p(x) / (x -a )

Tout se ramène donc au calcul d'une racine a d'un polynôme p(x)

2) Calcul d'une racine

  On cherche en partant de z   un nombre complexe a où la fonction  m(z)=|p(z)| admet un minimum local. D'après le second théorème  admis, a est alors une racine de p ( x ) .
Le calcul de a est résumé par l'algorithme suivant:

r = 1            # r est le pas initial
tant que r > 2 - 53 :
    si m(z) = 0, retourner z
    calculer une direction x = e it vers laquelle la fonction m(z) diminue        
    tant que m(z + r * x) < m(z):
        z = z + x
    r = r / 2
retourner z

  Le point délicat est le calcul d'une direction x = e i t vers laquelle la fonction m ( z ) diminue au voisinage de z .    
La formule de Taylor nous permet d'y arriver. Notons z = z 0 et k le plus petit entier ≥ 1 tel que P ( k ) ( z 0 ) 0 .
  Alors (Taylor) P ( z ) - P ( z 0 ) z 0 P ( k ) ( z 0 ) k ! ( z - z 0 ) k , donc P ( z ) z 0 P ( z 0 ) ( 1 + P ( k ) ( z 0 ) k ! P ( z 0 ) ( z - z 0 ) k ) .

  Ici, comme q ( x ) n'a que des racines simples, k = 1 , donc P ( z ) z 0 P ( z 0 ) ( 1 + P ' ( z 0 ) P ( z 0 ) ( z - z 0 ) ) .

  Donc m ( z ) m ( z 0 ) m ( 1 + P ' ( z 0 ) P ( z 0 ) ( z - z 0 ) ) . Si l'on part dans la direction e i t , alors z - z 0 = r e it (avec r > 0 ), donc   m ( z ) m ( z 0 ) m ( 1 + P ' ( z 0 ) P ( z 0 ) r e i t ) .

  Pour  que la fonction m diminue lorsque r croît à partir de r = 0 , il suffit que   P ' ( z 0 ) P ( z 0 ) r e i t soit un réel négatif, c'est à dire que   arg ( P ' ( z 0 ) P ( z 0 ) e i t ) = π [ 2 π ] , c'est à dire que t = - Arg ( P ' ( z 0 ) P ( z 0 ) ) + π .
  En pratique, on perturbe légèrement cette direction idéale, pour éviter des "impasses" dans certains cas particuliers.  On prendra donc t = - Arg ( P ' ( z 0 ) P ( z 0 ) ) + π + 0.01 et x = e i t .

Q3) Ecrire une fonction arg(z) (avec z nombre complexe non nul), calculant l'argument de z compris dans l'intervalle ]-π, π]
  On utilisera la fonction atan2(y, x) (dans le module math) qui calcule l'argument de z=x+iy compris dans l'intervalle ]-π, π].
  Par exemple: arg(-1)/pi, arg(1j)/pi, arg(-1 + sqrt(3)*1j)/pi = 1.0, 0.5, 0.6666666666666667

Q4) Ecrire une fonction direction(p, z) (avec p polynôme et z complexe) calculant le complexe x = e i t avec t = - arg ( P ' ( z ) P ( z ) ) + π + 0.01 .
  Par exemple, direction([1, 2, 3, 4, 5], 1 + 1j) = (-0.8349771271994405-0.5502846509341954j)

Q5) Pour rendre plus lisible les résultats, on décide d'arrondir les nombres complexe de la manière suivante:
  Ecrire une fonction arrondir(z) (avec z complexe) qui renvoie: (on note z = x + iy  et eps = 2**(-30)) :
    (1) 0 si |z| < eps
    (2) sinon x si |y| < eps
    (3) sinon i y si |x| < eps
    (4) sinon z

Par exemple, arrondir (-2. + 2.3e-17j) = -2.0 et arrondir(1j * z) = -2j .

Q6) Ecrire une fonction calcul_racine(p, z) (avec p polynôme) calculant une racine a de p en partant de z et en suivant les consignes données en III)2). On arrondira cette racine a avec la fonction arrondir(z).

  Par exemple, calcul_racine ( [ 1 , 1 , 1 ] , 3 + 1 j ) = ( - 0.5000000000000001 - 0.8660254037844387 j ) et calcul_racine ( [ 1 , 1 , 0 , 1 ] , 3 ) = - 0.6823278038280194

Q7) Ecrire une fonction liste_racines(p) (avec p polynôme) calculant la liste des racines du polynôme p en suivant les consignes données en III)1)

  Par exemple, liste_racines([1, 2, 3, 1j]) = [(-0.3945482804954703+0.43076327505918394j), (-0.27524656953737725-0.4752123804254154j), 0.6697948500328476+3.0444491053662315j)]

  Vérifier que cette fonction liste_racines(p) renvoie aussi  un résultat pour les polynômes avec des racines multiples, comme p = puissance ( [ 1 , 1 ] , 6 ) , mais de très mauvaise “qualité”.

IV) Phase (3) : calcul des racines du polynôme p(x) avec leur ordre de multiplicité

Q8) Ecrire une fonction ordre_multiplicite(p, z) qui calcule l'ordre n de la racine z du polynôme p. Si a n'est pas racine, n = 0. Utiliser le troisième théorème énoncé dans le I)
  Par exemple, ordre_multiplicite(p, -1/3) = 6 avec puissance([-5, -14, 3], 6).

Q9) Ecrire une fonction racines(p) (avec p polynôme) qui calcule la liste des racines de p répétées avec leur ordre de multiplicité.
Pour cela, calculer q = polynome_racines_simples(p), puis la liste s des racines de q et enfin exploiter la fonction ordre_multiplicite(p, z) pour terminer.

Par exemple:

Avec p = puissance([-5, -14, 3], 4), racines(p) = [-0.3333333333333333, -0.3333333333333333, -0.3333333333333333, -0.3333333333333333, 5.0, 5.0, 5.0, 5.0].

Tester avec p = produit ( puissance ( [ - 1 , 1 ] , 5 ) , puissance ( [ - 2 , 1 ] , 6 ) ) , avec p = produit ( puissance ( [ - pi , 1 ] , 5 ) , puissance ( [ - e , 1 ] , 6 ) ) et constater que les résultats sont excellents.

Tester avec p = puissance ( pol1 ( n ) , 2 ) pour n = 4, 5, 6, 7, 8 et constater la dégradation progressive de la qualité des résultats, voire même le flop pour n = 8 (que s'est-il passé ?)

V) Bilan

  Pour des polynômes de degrés raisonnables, avec des coefficients pas trop grands, et pas trop de racines multiples, les résultats sont excellents.
Une question importante et délicate serait de réfléchir aux points cruciaux des algorithmes écrits qui peuvent rendre les résultats tout à fait faux.

Q10) Quels sont ces points faibles ?

© 2014 - Eric Obermeyer      Powered by      Corrigé