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.

Pour un problème de recherche de racines, la fractale de Newton s'impose !

Les coefficients du polynôme de degré 10, bien qu'ils entrent dans la gamme des réels en double précision, conduisent à des dépassements de capacité lors des itérations. J'ai donc divisé tous les coefficients par 10000, ce qui ne change pas les racines.

Par ailleurs, une perturbation de 0,01% peut être appliquée au coefficient de x^9.

Voici d'abord le programme :

Code:: ' **********************************************************************
' Fractale de Newton : polynome de Wilkinson
' **********************************************************************

' Parametres de base (a modifier eventuellement)

const p = 10 ' Degre du polynome

data 3628800
data -10628640
data 12753576
data -8409500
data 3416930
data -902055
data 157773
data -18150
data 1320
data -55
data 1

const PicWidth = 600 ' Largeur de l'image
const PicHeight = 600 ' Hauteur de l'image

dim x0 = 5 ' Coord. X du centre
dim y0 = 0 ' Coord. Y du centre
dim MaxIter% = 200 ' Nb maxi d'iterations
dim ZoomFact = 0.2 ' Facteur de zoom
dim ColorFact = -2 ' Positif pour dessiner des bandes

const InsideCol = CL_NOIR ' Couleur des points non-convergents

' **********************************************************************

' Parametres supplementaires

const HalfPicWidth = PicWidth \ 2
const HalfPicHeight = PicHeight \ 2

dim ColFact, AbsCol, ScaleFact

SetParams()

' **********************************************************************

' Programme principal

dim Coef@(p), i%, coef_x

for i = 0 to p
  read coef_x
  coef(i) = cmplx(coef_x / 10000, 0)
next i

' Perturbation
' coef(9) = (1 + 0.0001) * coef(9)

dim x%, y%, btn%

mode 3, "Clic gauche pour zoomer / Clic droit pour dezoomer / ESC pour quitter", PicWidth, PicHeight

colormap PicWidth, PicHeight, fadr(newton)

repeat
  get_mouse x, y, btn

  if btn = 1 or btn = 2 then

   getcoord x0, y0, x, y

   SetZoom(btn = 1)
   SetParams()

   colormap PicWidth, PicHeight, fadr(newton)

  end_if

until inkey() = "ESCAPE"

' **********************************************************************

' Sous-programmes

sub init_sub (xt, yt, z@)
' Initialise les iterations au point (xt, yt)

  z = cmplx(xt, yt)
end_sub

sub iter_sub (z@, f@, df@)
' Calcul d'une iteration

  CPoly_Deriv z, coef(), f, df
end_sub

sub SetParams ()
  ColFact = 0.01 * ColorFact
  AbsCol = abs(ColFact)
  ScaleFact = 4 / (PicHeight * ZoomFact)
end_sub

sub SetZoom (zoom%)
  if zoom then
   ZoomFact = 5 * ZoomFact
   MaxIter = 1.25 * MaxIter
   ColorFact = 1.1 * ColorFact
  else
   ZoomFact = ZoomFact / 5
   MaxIter = MaxIter / 1.25
   ColorFact = ColorFact / 1.1
  end_if
end_sub

function rgbcol% (iter%, q, v)
' Determine la teinte (Hue) et la saturation d'apres le "Continuous Dwell"

  dim angle, radius, h, s, r%, g%, b%

  if q < 0.5 then
   q = 1 - 1.5 * q
   angle = 1 - q
  else
   q = 1.5 * q - 0.5
   angle = q
  end_if

  radius = sqr(q)

  if (ColFact > 0) and odd(iter) then
   v = 0.85 * v
   radius = 0.667 * radius
  end if

  h = frac(angle * 10) * 360
  s = frac(radius)

  HSVtoRGB h, s, v, r, g, b
  return RGB(r, g, b)
end_function

function newtcol% (iter%)
' Coloration pour la methode de Newton

  dim q, v

  v = 1
  q = log(iter) * AbsCol
  return rgbcol(iter, q, v)
end function

function newton% (x%, y%)
' Iteration de la fonction complexe au point (x, y)

  const Eps = 1.0E-6

  dim iter%
  dim z@, f@, df@, q@
  dim xt, yt

  xt = x0 + ScaleFact * (x - HalfPicWidth)
  yt = y0 + ScaleFact * (y - HalfPicHeight)

  init_sub xt, yt, z

  iter = 0

  repeat
   iter_sub z, f, df ' Calcul de f(z) et f'(z)

   if df = C_zero then return InsideCol

   q = f / df
   z = z - q ' z - f(z)/f'(z)

   iter = iter + 1
  until (cabs(f) < Eps and cabs(q) < Eps) or (iter > MaxIter)

  if iter > MaxIter then return InsideCol

  return newtcol(iter)
end_function

sub getcoord (x0, y0, x%, y%)
' Calcule les coordonnees (x0,y0) du centre
' d'apres la position (x,y) de la souris

  x0 = x0 + ScaleFact * (x - HalfPicWidth)
  y0 = y0 + ScaleFact * (y - HalfPicHeight)
end_sub

Et voici les résultats : en haut, sans perturbation, on retrouve bien nos 10 racines réelles sagement alignées sur l'axe horizontal ; en bas, avec perturbation il n'y a plus que 4 racines réelles (légèrement déplacées) ; les 6 autres ont éclaté en 3 groupes de 2 racines complexes conjuguées, situées symétriquement de part et d'autre de l'axe horizontal.

Avec l'ensemble de Mandelbrot c'est pareil : il faut diviser les coefficients sinon les "orbites" divergent très vite et on ne voit pas d'ensemble.

Un facteur autour de 10000 semble optimal. En-dessous l'ensemble se fragmente ; au-dessus il se ratatine et tend vers la figure correspondant à la fonction z^10 (voir l'exemple mandelmult.bas dans exemples\fractal).

Voici le programme :

Code:: ' ***************************************************
' Ensemble de Mandelbrot ou Julia : fonction polynome
' ***************************************************

' Parametres de base (a modifier eventuellement)

const PicWidth = 640 ' Largeur de l'image
const PicHeight = 480 ' Hauteur de l'image
const InsideCol = CL_BLANC ' Couleur de la partie interne

dim Nom$ ' Nom de l'image
dim p ' Degre du polynome
dim ICP% ' Indice du point critique
dim x0 ' Coord. X du centre
dim y0 ' Coord. Y du centre
dim MaxIter% ' Nb maxi d'iterations
dim ZoomFact ' Facteur de zoom
dim DistFact ' Plus negatif ==> contour plus sombre
dim ColorFact ' 0 pour noir et blanc, positif pour bandes de couleur
dim cx, cy ' Constante c pour l'ensemble de Julia
dim Julia& ' Ensemble de Julia ?

' -----------------------------------------------------
' Exemple : polynome de Wilkinson a l'ordre 10
' W(z) = (x - 1) * (x - 2) * ... * (x - 10) / K
' K ~ 10^4 produit des images plus interessantes
' La perturbation des coefficients deforme l'image et
' fait apparaitre des racines complexes dans la derivee
' -----------------------------------------------------

init_params "wilkinson",10,1,10,0,200,0.5,-2,-3,0,0

data 3628800
data -10628640
data 12753576
data -8409500
data 3416930
data -902055
data 157773
data -18150
data 1320
data -55
data 1

dim i%, a, coef@(10), z@(10)

for i = 0 to 10
  read a
  coef(i) = Cmplx(a / 10000, 0)
next i

' Perturbation
' coef(9) = (1 + 0.0001) * coef(9)

' ***************************************************

' Parametres supplementaires

const HalfPicWidth = PicWidth \ 2
const HalfPicHeight = PicHeight \ 2

const Esc = 1.0E+10 ' Rayon d'echappement
const LnEsc = log(Esc)

dim Lnp : Lnp = log(p)

dim ColFact, AbsCol, ScaleFact

SetParams()

' ***************************************************

' Calcul des points critiques (racines de la derivee)

dim coefd@(p - 1), root@(p - 1)

coefd(0) = coef(1)
for i = 2 to p
  coefd(i - 1) = i * coef(i)
next i

CPoly_Roots coefd(), root()

? : ? "Points critiques"
? "---------------------------"
for i = 1 to p - 1
  ? using "z(#) = "; i;
  ? using "####.####"; root(i).x;
  if root(i).y <> 0 then ? using "+####.####i"; root(i).y;
  ?
next i

' Dessin de l'image

dim x%, y%, btn%, key$

mode 3, "Clic gauche : zoom + / Clic droit : zoom - / S : sauvegarde / ESC : quitte", PicWidth, PicHeight

colormap PicWidth, PicHeight, fadr(mandelbrot)

repeat
  get_mouse x, y, btn

  if btn = 1 or btn = 2 then
   getcoord x0, y0, x, y

   SetZoom(btn = 1)
   SetParams()

   colormap PicWidth, PicHeight, fadr(mandelbrot)
  end_if

  key = inkey()
  if upper(key) = "S" then save()
until key = "ESCAPE"

' ***************************************************

' Sous-programmes

sub init_params(_Nom$, _p, _ICP%, _x0, _y0, _MaxIter%, _ZoomFact, _DistFact, _ColorFact, _cx, _cy)
' Initialise les parametres

  Nom = _Nom
  p = _p
  ICP = _ICP
  x0 = _x0
  y0 = _y0
  MaxIter = _MaxIter
  ZoomFact = _ZoomFact
  DistFact = _DistFact
  ColorFact = _ColorFact
  cx = _cx
  cy = _cy

  Julia = (cx <> 0) or (cy <> 0)
end_sub

sub init_sub (xt, yt, c@, z@, dz@)
' Initialise les iterations au point (xt, yt)

  if Julia then
   c = cmplx(cx, cy)
   z = cmplx(xt, yt)
   dz = C_one
  else
   c = cmplx(xt, yt)
   z = root(ICP)
   dz = C_zero
  end if
end_sub

sub iter_sub (c@, z@, dz@, zn@, dzn@)
' Calcul d'une iteration

  dim f@, df@

  CPoly_Deriv z, coef(), f, df

  zn = f + c ' f(z) + c
  dzn = df * dz ' Derivee : dz/dc

  if not Julia then dzn = dzn + 1
end_sub

sub SetParams ()
  ColFact = 0.01 * ColorFact
  AbsCol = abs(ColFact)
  ScaleFact = 4 / (PicHeight * ZoomFact)
end_sub

sub SetZoom (zoom%)
  if zoom then
   ZoomFact = 5 * ZoomFact
   MaxIter = 1.25 * MaxIter
   ColorFact = 1.1 * ColorFact
  else
   ZoomFact = ZoomFact / 5
   MaxIter = MaxIter / 1.25
   ColorFact = ColorFact / 1.1
  end_if
end_sub

function rgbcol% (iter%, q, v)
' Determine la teinte (Hue) et la saturation d'apres le "Continuous Dwell"

  dim angle, radius, h, s, r%, g%, b%

  if q < 0.5 then
   q = 1 - 1.5 * q
   angle = 1 - q
  else
   q = 1.5 * q - 0.5
   angle = q
  end_if

  radius = sqr(q)

  if (ColFact > 0) and odd(iter) then
   v = 0.85 * v
   radius = 0.667 * radius
  end if

  h = frac(angle * 10) * 360
  s = frac(radius)

  HSVtoRGB h, s, v, r, g, b
  return RGB(r, g, b)
end_function

function mdbcol% (iter%, mz, mdz)
' Coloration pour les ensembles de Mandelbrot/Julia
' Retourne la couleur RGB d'un point en fonction de :
' iter = nb d'iterations
' mz, mdz = modules de z et de (dz/dc) a l'iteration iter
' Methode de coloration d'apres R. Munafo (http://mrob.com/pub/muency/color.html)

  dim lmz, dist, dwell, dscale, q, v

  ' Determiner la luminosite (Value) d'apres l'estimateur de distance

  lmz = log(mz)

  dist = p * mz * lmz / mdz
  dscale = log(dist / ScaleFact) / Lnp + DistFact

  if dscale > 0 then
   v = 1
  elseif dscale > -8 then
   v = 1 + dscale / 8
  else
   v = 0
  end_if

  dwell = iter + log(LnEsc / lmz) / Lnp
  q = log(abs(dwell)) * AbsCol
  return rgbcol(iter, q, v)
end_function

function mandelbrot% (x%, y%)
' Iteration de la fonction complexe au point (x, y)

  dim iter%
  dim c@, z@, dz@, zn@, dzn@
  dim xt, yt, mz, mdz

  xt = x0 + ScaleFact * (x - HalfPicWidth)
  yt = y0 + ScaleFact * (y - HalfPicHeight)

  init_sub xt, yt, c, z, dz

  iter = 0
  mz = cabs(z)

  while iter < MaxIter and mz < Esc
   iter_sub c, z, dz, zn, dzn
   z = zn
   dz = dzn
   mz = cabs(z)
   iter = iter + 1
  wend

  if iter = MaxIter then return InsideCol

  mdz = cabs(dz)
  return mdbcol(iter, mz, mdz)
end_function

sub getcoord (x0, y0, x%, y%)
' Calcule les coordonnees (x0,y0) du centre
' d'apres la position (x,y) de la souris

  x0 = x0 + ScaleFact * (x - HalfPicWidth)
  y0 = y0 + ScaleFact * (y - HalfPicHeight)
end_sub

sub save ()
' Sauvegarde l'image et les parametres

  img_save Nom + ".png"
  openout #1, Nom + ".pol"
  write #1, Nom, p, ICP, x0, y0, MaxIter, ZoomFact, DistFact, ColorFact, cx, cy
  closeout #1
end_sub

Et voici le résultat : en haut, l'ensemble non perturbé, en bas : l'ensemble perturbé, avec à chaque fois comme point critique la première racine de la dérivée qui est dans les deux cas une racine réelle.

Une perturbation de 0,01% sur le coefficient de z^9 suffit à fragmenter complétement l'ensemble !

» 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