Recherche :

Calculs des zéros de polynômes
Les réponses
Licence de mathématiques Nouméa 2004

Jean-Philippe CHANCELIER
(avec le soutien du Ministère de l'Outre-Mer)
Version pdf de ce document
Version sans bandeaux

Table des matières

1 Calcul des valeurs d'un polynôme et de ses dérivées par l'algorithme de Horner

Soit $ P(x)$ un polynôme de degré $ n$ :

$\displaystyle P(x) \stackrel{\mbox{\tiny def}}{=}\sum_{i=0}^n a_{i+1} x^i $

On cherche ici à évaluer la valeur du polynôme pour une valeur donnée $ y$. Pour ce faire, on peut factoriser le polynôme $ P(x)$ sous la forme :

$\displaystyle P(x) = Q(x)(x-y) + R$ (1)

où $ Q(x)$ est un polynôme de degré $ n-1$ et $ R$ une constante scalaire. Sous cette forme on voit que $ P(y)=R$ et que $ P'(y)=Q(y)$. Le polynôme $ Q(x)$ est bien sur différent du polynôme $ P'(x)$ mais leur valeurs coïncident quand $ x$ prends la valeur $ y$.

Soit $ b_i$ les coefficients du polynôme $ Q$ :

$\displaystyle Q(x) \stackrel{\mbox{\tiny def}}{=}\sum_{i=0}^{n-1} b_{i+1} x^i$ (2)

en développant l'équation (1) on obtient facilement une équation récurrente permettant de calculer les coefficients du polynôme $ Q$ :

$\displaystyle b_i = b_{i+1}y + a_{i+1},\quad b_{n}= a_{n+1},$   pour$\displaystyle \quad i= n-1,n-2,\ldots, 0$ (3)

La valeur de $ R$ est donnée par $ b_0$. On notera que l'on peut initialiser la récurrence par $ b_{n}= a_{n+1}$ ou quitte à rajouter un point par $ b_{n+1}=0$.

Cet Algorithme s'appelle algorithme d'Horner, il revient à évaluer la valeur du polynôme $ P$ en le factorisant sous la forme :

$\displaystyle P(x) = ((\cdots (a_{n+1}*x + a_{n})*x + a_{n-1})*x + \cdots)*x + a_1 $

Nous avons noté que la dérivée de $ P(x)$ au point $ y$ s'obtient en évaluant $ Q(y)$. Il suffit donc de réappliquer l'algorithme d'Horner au polynôme $ Q$ pour évaluer $ P'(y)$.


->function [Q,R]=Horner(P,xi)   $ \leftarrow$ Évalue le polynôme P au point xi
->  [lhs,rhs]=argn(0);
->  a=coeff(P);   $ \leftarrow$ Les coefficients du polynôme
->  na=size(a,'*');
->  b=ones(a);    $ \leftarrow$ Espace pour stocker les coefficient du polynôme Q et la valeur R=P(xi)
->  b($)=a($);    $ \leftarrow$ Initialisation
->  for i=na-1:-1:1
->    b(i)= b(i+1)*xi+ a(i)   $ \leftarrow$ Les coefficients du polynôme Q
->  end
->  if lhs==2 then
->    R = b(1);
->    Q=poly(b(2:$),'x','coeffs');   $ \leftarrow$ Construction de Q
->  else
->    Q= b(1)
->  end
->endfunction

->P=poly(1:3,'x','coeffs');   $ \leftarrow$ $ P(x)= 1 + 2x + 3 x\sp{2} $

->xi=0.5;
->[Q,R]=Horner(P,xi);

->P2= Q*poly(xi,'x','roots')+R   $ \leftarrow$ On vérifie le lien entre P, Q et R
 P2  =

               2  
    1 + 2x + 3x   

->if norm(coeff(P-P2)) > %eps then pause;end

->if norm( R - horner(P,xi)) > %eps then pause;end   $ \leftarrow$ Contrôle avec la fonction Scilab horner

->[Pprimxi]=Horner(Q,xi);   $ \leftarrow$ Calcul de $ P'(xi)$

->if norm( Pprimxi - horner(derivat(P),xi)) > %eps then pause;end   $ \leftarrow$ Contrôle avec les fonctions Scilab horner et derivat

En fait, on peut calculer directement les valeurs $ P(y)$ et $ P'(y)$ sans passer par l'intermédiaire de la construction du polynôme $ Q$, c'est à dire sans stocker les coefficients du polynôme $ Q$. On écrit conjointement les deux algorithmes de Horner et pour faciliter la gestion des indices on initialise le deuxième avec $ c_n=0$ plutôt que $ c_{n-1}=b_n$. Ainsi les deux vecteurs de coefficients à calculer auront même taille :

$\displaystyle b_i$ $\displaystyle =$ $\displaystyle b_{i+1}y + a_{i+1},\quad b_{n}= a_{n+1},$   pour$\displaystyle \quad i= n-1,n-2,\ldots, 0$  
$\displaystyle c_i$ $\displaystyle =$ $\displaystyle c_{i+1}y + b_{i+1},\quad c_{n}= 0,$   pour$\displaystyle \quad i= n-1,n-2,\ldots, 0$  

le calcul du couple $ (b_0,c_0)$ permet d'obtenir le couple de valeurs $ (P(y),P'(y)$. En généralisant cette idée on peut calculer les $ p$-dérivées successives du polynôme $ P$ en un point $ y$. Posons $ P_i(x)= P_{i+1}(x)(x-y) + R_i$ avec $ P_0(x)=P(x)$. Par un algorithme de Horner on peut calculer en procédant comme dans (3) les valeurs $ P_i(y)$ pour $ i=0,\ldots,p-1$. Il est ensuite facile de voir par dérivation composées que :

$\displaystyle P^{(i)}_0(y) = i! P_{i}(y) ,$   pour$\displaystyle \quad i\ge 1$ (4)

et d'obtenir ainsi les dérivées successives de $ P$ en $ y$.


->function [val]=Hornerp(P,xi,p)   $ \leftarrow$ Évalue $ P\sp{(i)}(xi)$ pour $ i=0,\ldots,p-1)$
->  V=zeros(p,1);
->  a=coeff(P);
->  V(1)=a($);   $ \leftarrow$ Initialisation
->  M=diag(xi*ones(1,p))+ diag(ones(1,p-1),-1);
->  n=size(a,'*');
->  for i=1:n-1
->    V=M*V + [a(n-i);zeros(p-1,1)],   $ \leftarrow$ Horner simultané
->  end
->  F=ones(p,1);
->  for i=2:p-1 do F(i+1)=F(i)*i ; end   $ \leftarrow$ Calcul des $ i!$
->  val = V.* F;   $ \leftarrow$ Les dérivées de P en $ y$
->endfunction

->P=poly(1:8,'x','coeffs');
->xi=1.5 ;
->p=10;
->val=Hornerp(P,xi,p);

->Q=P;
->for i=1:p-1, Q=[Q;derivat(Q($))];end   $ \leftarrow$ Les dérivées succéssives de P

->if norm(val -horner(Q,xi)) > %eps then pause;end    $ \leftarrow$ Vérifions le résultat

1.1 Utilisation de l'algorithme de Horner pour la factorisation

Supposons que l'on ait estimé un zéro $ \xi$ du polynôme $ P$, dans la décomposition de $ P(x)=Q(x)(x-\xi)+R$ on doit avoir $ R=0$ et la méthode de Horner permet d'éliminer la racine $ \xi$ du polynôme $ P(x)$ puisqu'elle donne les coefficients du polynôme $ Q$. Mais si $ \xi$ est une racine de $ P(x)$ on doit aussi avoir $ R=b_0=0$ dans l'algorithme de Horner.

On peut donc utiliser la récurrence (3) à partir de de $ b_{n}= a_{n+1}$ (méthode forward) ou à partir de $ b_0=0$ méthode backward. Si l'on cherche à éliminer successivement les racines d'un polynôme en éliminant à chaque fois la plus grande racine, la méthode backward donne une meilleur stabilité numérique.

On peut comparer les deux méthodes avec le polynôme $ P(x)\stackrel{\mbox{\tiny def}}{=}\prod_{j=0}^{13}(x - 2^{-j})$ en utilisant la fonction roots de Scilab pour estimer les valeurs des racines.


->function [Q]=deflate_forward(P,xi)   $ \leftarrow$ Élimine la racine xi méthode forward
->  a=coeff(P);   $ \leftarrow$ Les coefficients du polynôme
->  na=size(a,'*');
->  b=ones(a);    $ \leftarrow$ Espace pour stocker les coefficient du polynôme Q et la valeur R=P(xi)
->  b($)=a($);    $ \leftarrow$ Initialisation
->  for i=na-1:-1:1
->    b(i)= b(i+1)*xi+ a(i)   $ \leftarrow$ Les coefficients du polynôme Q
->  end
->  Q=poly(b(2:$),'x','coeffs');   $ \leftarrow$ Construction de Q
->endfunction

->function [Q]=deflate_backward(P,xi)   $ \leftarrow$ Élimine la racine xi méthode backward
->  a=coeff(P);   $ \leftarrow$ Les coefficients du polynôme
->  na=size(a,'*');
->  b=ones(a);    $ \leftarrow$ Espace pour stocker les coefficient du polynôme Q et la valeur R=P(xi)
->  b(1)=0;    $ \leftarrow$ Initialisation
->  for i=1:na-1
->    b(i+1)= (b(i)-a(i))/xi $ \leftarrow$ Les coefficients du polynôme Q
->  end
->  Q=poly(b(2:$),'x','coeffs');   $ \leftarrow$ Construction de Q
->endfunction


->n=13;
->r= 2.^(-[0:n]);
->P=poly(r,'x','roots');
->rr=roots(P)
 rr  =

!   0.0001221 !
!   0.0002441 !
!   0.0004883 !
!   0.0009766 !
!   0.0019531 !
!   0.0039062 !
!   0.0078125 !
!   0.015625  !
!   0.03125   !
!   0.0625    !
!   0.125     !
!   0.25      !
!   0.5       !
!   1.        !

->maxi(abs(rr'-r($:-1:1)))
 ans  =

    4.441D-16  

->P1=P;
->r3=zeros(r);
->for i=1:n+1
->  r3(i)=maxi(real(roots(P1)));
->  P1=deflate_backward(P1,r3(i));    $ \leftarrow$ On élimine la plus grande racine
->end
->maxi(abs(r3-r))      $ \leftarrow$ Test du résultat
 ans  =

    6.661D-16  

->P1=P;
->r4=zeros(r);
->for i=1:n+1
->  r4(i)=maxi(real(roots(P1)));
->  P1=deflate_forward(P1,r4(i));   $ \leftarrow$ On élimine la plus grande racine
->end

->maxi(abs(r4-r))   $ \leftarrow$ Test du résultat
 ans  =

    0.0408097  

1.2 La méthode de Horner backward

L'équation (3) permettant de calculer les coefficients du polynôme $ Q$ peut-être vu comme un système dynamique discret. La stabilité de ce système ou de celui obtenu pour le calcul simultané de $ P$ et de ses dérivées au point $ y$ dépend de $ \vert y\vert$. Pour $ \vert y\vert>1$ le système est instable et on va donc pour des polynômes de grande taille amplifier les erreurs de calcul si l'on utilise (3). Mais on peut remarquer que $ P(x)$ s'écrit aussi :

$\displaystyle P(x) = x^n \tilde{P}(1/x),$   avec$\displaystyle \quad \tilde{P}(x)\stackrel{\mbox{\tiny def}}{=}\sum_{i=0}^n a_{i+1} x^{(n-i)} )$ (5)

Pour une valeur de $ y$ telle que $ \vert y\vert>1$. on se ramène à un système dynamique stable en utilisant l'algorithme de Horner pour l'évaluation de $ \tilde{P}$ et de ses dérivées en $ 1/y$ et on utilise la relation (5) pour les relier aux valeurs de $ P$ et de ses dérivées en $ y$.


->n=13;
->r= 2.^(-[0:n]);
->P=poly(r,'x','roots');

->xi=10000 ;
->p=2;
->val=Hornerp(P,xi,p);    $ \leftarrow$ $ P(xi)$ par une méthode forward

->xii=(1/xi);
->a=coeff(P);
->Ptilde=poly(a($:-1:1),'x','coeffs');   $ \leftarrow$ Construction de $ \tilde{P}$
->tval=Hornerp(Ptilde,xii,p);    $ \leftarrow$ $ \tilde{P}(1/xi)$ par une méthode forward

->n = degree(P);
->back(1)= tval(1)*xi^n;    $ \leftarrow$ $ P(xi)$
->back(2)= n*xi^(n-1)*tval(1) - xi^(n-2)*tval(2);    $ \leftarrow$ $ P'(xi)$

->norm((back-val)./val)
 ans  =

    4.614D-16  

2 Recherche des zéros d'un polynôme

On utilise ici la méthode de Newton pour trouver les zéros d'un polynôme. Soit $ P(x)$ un polynôme de degré $ n$, on utilise pour trouver les zéros l'algorithme de Newton qui consiste à effectuer les itérations suivantes :

$\displaystyle x_{k+1} = x_k - P(x_k)/P'(x_k)$ (6)

Il faut pouvoir évaluer la valeur du polynôme et de sa dérivée au point courant $ x_k$ et nous avons vu qu'un algorithme d'Horner permet d'évaluer ces deux quantités.

La convergence de l'algorithme de Newton dans ce cas particulier est donnée par le théorème suivant :

Théorème 1   Soit $ P$ un polynôme de degré $ n$ dont toutes les racines sont réelles. La méthode de Newton donne une suite strictement décroissante pour $ x_0 > \xi_1$ où $ xi_1$ est la plus grande racine de $ P$.

On trouvera la preuve de la convergence dans [1]. Pour implémenter cet algorithme il faut trouver un majorant de $ xi_1$. Plusieurs formules de majoration sont données dans la littérature. Par exemple :

$\displaystyle \vert\xi_1\vert < \max\left\{ 1, \sum_{j=0}^{n-1} \frac{\vert a_{j+1}\vert}{\vert a_{n+1}\vert}\right\}$   si$\displaystyle \quad P(x) = \sum_{i=0}^n a_{i+1} x^i$ (7)

$\displaystyle \vert\xi_1\vert < \max \left\{ \frac{\vert a_1\vert}{\vert a_{n+1... ...t a_{n+1}\vert},\ldots, 1+ \frac{\vert a_{n}\vert}{\vert a_{n+1}\vert} \right\}$   si$\displaystyle \quad P(x) = \sum_{i=0}^n a_{i+1} x^i$ (8)


->function xn=newton(P,x,prec)
->  while %t
->    val=Hornerp(P,x,2);    $ \leftarrow$ val=[P(x),P'(x)]
->    xn= x -val(1)/val(2);
->    if norm(xn-x) <= prec*norm(x) then
->      break;
->    end
->    x=xn;
->  end
->endfunction

->n=13;
->r= 2.^(-[0:n]);
->P=poly(r,'x','roots');    $ \leftarrow$ Le polynôme test

->a=coeff(P);
->xd=max([1+abs(a(2:$-1)/a($)),abs(a(1)/a($))]) ;   $ \leftarrow$ Un majorant de la plus grande racine

->y1=newton(P,xd,%eps)
 y1  =

    1.  

En éliminant à chaque fois la plus grande racine trouvée on peut chercher toutes les racines du polynôme $ P(x)$. On peut constater dans le code qui suit que la méthode d'élimination choisie (forward ou backward) n'est pas anodine !


 y1  =

    1.  

->y=y1; rnf= zeros(r); rnf(1)=y; Q=P;
->for i=2:4    $ \leftarrow$ Problème au dela de $ i=4$
->  Q=deflate_forward(Q,y); $ \leftarrow$ On élimine la plus grande racine
->  y=newton(Q,y,%eps);   $ \leftarrow$ Évaluation de la plus grande racine
->  rnf(i)=y;
->end

->y=y1;rnb= zeros(r); rnb(1)=y; Q=P;
->for i=2:(n+1)
->  Q=deflate_backward(Q,y);   $ \leftarrow$ On élimine la plus grande racine
->  y=newton(Q,y,%eps);
->  rnb(i)=y;
->end

->if norm(abs(r-rnb)) > 10*%eps then pause;end

Plutôt que d'essayer de simplifier le polynôme $ P$ on peut utiliser la méthode de Maehly (1954) qui consiste à appliquer la méthode de Newton à la fraction rationnelle $ P(x)/(x-\xi)$ plutôt que d'éliminer la racine $ \xi$ du polynôme $ P(x)$. Ainsi si l'on a estimé les $ j$ premières racines du polynôme, notées $ \xi_i$, on cherche la $ j+1$-ème racine par l'algorithme de Newton :

$\displaystyle x_{k+1}= x_k - Q(x_k)$   avec$\displaystyle \quad Q(x) \stackrel{\mbox{\tiny def}}{=}\frac{P(x)}{ P^{(1)}(x) - \sum_{i=1}^{j} \frac{P(x)}{x - \xi_i}}$ (9)

Pour initialiser l'algorithme, on peut utiliser la dernière racine trouvée (à epsilon près pour éviter une division par zéro).


->function y=newton(P,x,r,prec)
->  while %t
->    val=Hornerp(P,x,2);
->    if r<>[] then
->      w=sum(val(1)./(x-r)); $ \leftarrow$ r $ =\sum\sb{i=1}\sp{j} \frac{P(x)}{x - \xi\sb{i}}$
->    else
->      w=0.0;
->    end
->    xn= x -val(1)/(val(2)-w);
->    if norm(xn-x) <= prec*norm(x) then
->      break;
->    end
->    x=xn;
->  end
->  y=xn;
->endfunction

->n=13;r= 2.^(-[0:n]);P=poly(r,'x','roots');
->y=newton(P,10,[],%eps);

->rn= zeros(r); rn(1)=y;

->for i=2:n+1
->  y=newton(P,y-100*%eps,rn(1:i-1),10*%eps);
->  rn(i)=y;
->end

->if maxi(abs(rn-r)) > 10*%eps then pause;end

Bibliographie

1
J. Stoer and R. Burlisch.
Introduction to Numerical Analysis.
Springer-Verlag, 1980.