Le polynôme de Wilkinson

La dernière version de FBCroco inclut une fonction CPoly_Roots pour calculer les racines réelles ou complexes d'un polynôme à coefficients réels ou complexes.

Le polynôme de Wilkinson constitue un exemple classique pour tester de tels programmes. C'est un polynôme du vingtième degré formé par un produit de facteurs (x - i) où i est un entier de 1 à 20.

Code:: W(x) = (x - 1) * (x - 2) * ... * (x - 20)

Les nombres de 1 à 20 sont donc les racines du polynôme.

En développant W(x) sous la forme a0 + a1 * x + a2 * x^2 + ... on constate que les coefficients peuvent avoir jusqu'à 20 chiffres (voir les valeurs dans les DATA du programme ci-dessous). C'est plus que ce que peuvent traiter les réels "double précision" qui sont limités à environ 15 chiffres. On va donc passer en précision étendue avec MPFR, en prenant 30 chiffres pour avoir une marge.

Code:: #set_mpfr_prec 30

Les coefficients sont lus sous forme de chaînes de caractères, automatiquement traduites en nombres complexes par la fonction Cmplx. On appelle ensuite la fonction CPoly_Roots. Pour avoir les racines avec 10 chiffres exacts on va régler la tolérance à 10^(-12) et le nombre maximal d'itérations à 100, d'où l'instruction :

Code:: CPoly_Roots coef(), z(), 100, 1E-12

Voici le programme complet :

Code:: ' *******************************************************
' Calcul des racines du polynome de Wilkinson
' P(x) = (x-1)(x-2)...(x-20)
' https://fr.wikipedia.org/wiki/Polynôme_de_Wilkinson
' *******************************************************

#set_mpfr_prec 30

data " 2432902008176640000"
data " -8752948036761600000"
data " 13803759753640704000"
data " -12870931245150988800"
data " 8037811822645051776"
data " -3599979517947607200"
data " 1206647803780373360"
data " -311333643161390640"
data " 63030812099294896"
data " -10142299865511450"
data " 1307535010540395"
data " -135585182899530"
data " 11310276995381"
data " -756111184500"
data " 40171771630"
data " -1672280820"
data " 53327946"
data " -1256850"
data " 20615"
data " -210"
data " 1"

dim i%, a$, coef@(20), z@(20)

for i = 0 to 20
  read a
  Coef(i) = Cmplx(a, 0)
next i

CPoly_Roots Coef(), z(), 100, 1E-12

cls
? "Racines du polynome de Wilkinson"
? "--------------------------------"

for i = 1 to 20
  ? using "##.########"; mpfr_to_dbl(z(i).x)
next i

Et voici le résultat (tel qu'il apparaît, dans le désordre) :

Code:: Racines du polynome de Wilkinson
--------------------------------
1.00000000
8.00000000
16.00000000
19.00000000
11.00000000
5.00000000
6.00000000
14.00000000
20.00000000
13.00000000
7.00000000
2.00000000
10.00000000
18.00000000
17.00000000
9.00000000
3.00000000
4.00000000
12.00000000
15.00000000

Salut tout le monde.
Salut Jean.

En voyant ce sujet, quelques souvenirs se sont bousculés dans ma tête.
J'ai déjà vu quelque chose de ce genre, mais où et quand ?
Heureusement la maladie qui m'accompagne depuis un certain temps n'a pas pu influencer ma mémoire.

Bon, j'ai retrouvé : Octobre 1985 - SOFT & MICRO
Je cite le magazine :

"Soit P le polynôme suivant :
P(x) = (x+1).(x+2)...(x+20) = x²⁰ + 210.x¹⁹ + ... + 20!
C'est un polynôme de degré 20, et ses racines sont évidemment les 20 entiers :
-1, -2, ..., -20
Ajoutons maintenant au coefficient de x¹⁹ une très faible perturbation : 2^-23.
On pourrait s'attendre à ce que la disposition des racines ne change que très peu.
Le résultat est surprenant."
Fin de citation.

On voit que les 9 premières solutions et la 20ème restent les mêmes (inchangées) alors que les dix autres solutions (de 10 à 19) sont devenues des racines complexes conjuguées deux à deux.

Je cite le magazine:
" L'analyse de ce phénomène relève de la théorie des catastrophes : celle-ci a été récemment édifiée par le mathématicien français René Thom pour rendre compte de systèmes tels qu'une faible perturbation de certains de leurs paramètres entraine une violente variation qualitative de leur comportement.
Fin de citation.

Merci papydall. C'est exactement ça !

Pour mieux comprendre, traçons le graphe du polynôme selon la représentation semi-logarithmique utilisée sur la page de Wikipedia :

Code:: ' *******************************************************
' Calcul des racines du polynome de Wilkinson
' W(x) = (x-1)(x-2)...(x-20)
' https://fr.wikipedia.org/wiki/Polynôme_de_Wilkinson
' *******************************************************

#set_mpfr_prec 30

data " 2432902008176640000"
data " -8752948036761600000"
data " 13803759753640704000"
data " -12870931245150988800"
data " 8037811822645051776"
data " -3599979517947607200"
data " 1206647803780373360"
data " -311333643161390640"
data " 63030812099294896"
data " -10142299865511450"
data " 1307535010540395"
data " -135585182899530"
data " 11310276995381"
data " -756111184500"
data " 40171771630"
data " -1672280820"
data " 53327946"
data " -1256850"
data " 20615"
data " -210"
data " 1"

dim i%, a$, coef@(20), z@(20)

for i = 0 to 20
  read a
  Coef(i) = Cmplx(a, 0)
next i

' Perturbation de 2^(-23) sur le coefficient de x^19
' Coef(19) = Coef(19) - mpfr(2)^(-23)

' Calcul des racines

CPoly_Roots Coef(), z(), 100, 1E-12

cls
print "Racines du polynome de Wilkinson"
print "--------------------------------"
print " Real Imag. "

for i = 1 to 20
  print using "####.########"; z(i).x;
  if abs(z(i).y) > 1E-10 then
   print using "####.########"; z(i).y
  else
   print
  end if
next i

' Graphe du polynome (cf Wikipedia)

mode 2, "Polynome de Wilkinson"

fplot fadr(func), "x", "LIN 0 20 2", "sgn(W) * ln(1 + abs(W))", "LIN -40 40 10", "0 20 18"

while inkey = "" : wend

end

function func(x)
  dim z@, p@
  z = cmplx(x, 0)
  p = CPoly(z, Coef())
  return sgn(p.x) * log(1 + abs(p.x))
end_function

Le graphe du polynôme non perturbé :

L'échelle verticale étant logarithmique (log népérien), on voit qu'au voisinage des racines le polynôme passe brutalement de e^(-40) à e^40,
soit environ 34 ordres de grandeurs !

Pour le polynôme perturbé, on voit la disparition d'une partie des racines réelles :

Pour le polynôme perturbé les racines calculées sont les suivantes :

Code:: Racines du polynome de Wilkinson
--------------------------------
   Real Imag.
   1.00000000
   8.00726760
  16.73073747 2.81262489
  19.50243940 -1.94033035
  10.09526615 -0.64350090
   4.99999993
   6.00000694
  13.99235814 2.51883007
  20.84690810
  11.79363388 -1.65232973
   6.99969723
   2.00000000
  10.09526615 0.64350090
  19.50243940 1.94033035
  16.73073747 -2.81262489
   8.91725025
   3.00000000
   4.00000000
  11.79363388 1.65232973
  13.99235814 -2.51883007

Autre essai avec le polynôme de degré 10. Ici on n'a plus besoin de la précision étendue.

La perturbation sur le coefficient de x^9 (soit -55) est de 2^(-10) soit environ 0.001.

Les résultats sont analogues à ceux obtenus avec le polynôme de degré 20.

Il sera intéressant de voir ce que cela donne en termes d'images fractales.

» Evaluer un polynôme en un point par la méthode de Horner
» Calcul des zéros réels et / ou complexes d’un polynôme