Électro-localisation de racines

i.zitoussi

Salut

Un petit fil pour se détendre, qui intéressera peut-être ceux qui ne connaissent pas déjà.

J'ai essayé pour la n-ième fois de faire un algorithme de factorisation de polynômes, et c'est vraiment dur. Il faut trouver une racine, factoriser, recommencer, et les erreurs d'approximations peuvent s'amplifier rapidement.

Et puis j'ai lu cette page wikipedia sur un procédé vraiment magique et simple comme bonjour, qui peut se résumer en : pour factoriser (numériquement, c'est-à-dire en admettant des approximations) un polynôme unitaire complexe $P$ de degré $N$, il suffit de placer $N$ électrons dans le plan complexe muni d'un potentiel électrique égal à $\frac{P'}{P}$, et puis ... attendre. Pas très longtemps, car ça marche remarquablement bien.

Je ne donne pas les détails car la page wikipedia le fait assez bien, et l'article original de l'auteur de l'algorithme (O. Aberth) cité dans cette page est simple et en libre accès.

J'ai fait quelques essais en prenant tout au hasard dans le carré $|\Re(z)|\leq 5$, $|\Im(z)|\leq 5$ :

• les cercles rouges sont les protons/racines du polynôme (qui ne bougent pas) ;
• les cercles bleus sont les électrons/valeurs d'essai initiales des racines (appelées $z_i$, $i=1,\ldots,N$, dans la suite) ;
• les trajectoires turquoises sont celles des "électrons" fournies par l'algorithme de Newton:$$z_i\mapsto z_i-\frac{P(z_i)}{P'(z_i)} \qquad (\mathrm{Newton})$$ c'est-à-dire "d'un électron comme s'il était seul dans ce champ" ;
• et les trajectoires bleues sont celles fournies par l'algorithme:$$z_i\mapsto z_i-\frac{1}{\frac{P'(z_i)}{P(z_i)} + V_i(z_i)}\qquad (\mathrm{Aberth})$$ où $V_i(z_i)$ est le potentiel dû à la présence des "autres électrons"...

Le résultat est étonnant.

Dans l'algorithme de Newton, le procédé est régulier mais relativement lent, et plusieurs "électrons" peuvent choisir le même proton.

L'algorithme de Aberth (MPSolve) est plus chaotique, mais nécessite beaucoup moins d'itérations, et chaque électron se trouve son propre proton. Si deux électrons se dirigent vers le même proton, fatalement l'un d'entre eux se fait éjecter, et quitte à sortir de la scène il revient presque aussitôt pour se trouver un autre proton libre.

Une capture d'écran pour un polynôme aléatoire de degré 100 (les deux algos sur une même figure), et en attachement, les deux algos séparés:

raoul.S

Vraiment géniale l'idée. Je me suis amusé à suivre quelques trajectoires sur tes images et effectivement l'algorithme d'Aberth est vraiment mieux pour chopper toutes les racines.

PS. il aurait peut-être fallu poster dans Analyse, car on est loin du Shtam quand-même

Calli

Bonjour,
Merci pour le partage !

i.zitoussi

Raoul, c'est le Shtam utilisé comme section "Maths Générales", vu qu'il n'y en a pas. Mais j'ai mis ça là aussi parce qu'il y a quelques jours il y a eu ce fil  dans le Shtam sur les zéros de la fonction de Riemann versus "certains phénomènes physiques" (comprendra qui le peut, pas moi en tout cas), et ça m'y a fait penser. J'avais d'ailleurs oublié d'y faire référence dans mon post...

Sinon, on est bien d'accort, c'est génial. Wikipedia dit que l'analogie $\frac{P'}{P} \leftrightarrows$ potentiel de Coulomb dû à $N$ charges égales revient à Stieltjes, mais ne dit pas ce qu'il en a fait. L'algorithme, lui, date de 1969-1973 (deux auteurs en fait, et séparés, pas trop compris qui a fait quoi).

Le point vraiment super est qu'une fois "qu'un électron s'est approprié un proton", c-à-d que sa position tend à se confondre avec celle d'un proton-racine, la couple proton/électron disparaît des radars pour les autres électrons (le champ résultant est nul), et c'est comme si la factorisation se faisait toute seule, sans perte de précision !

Calli

i.zitoussi a dit :

Le point vraiment super est qu'une fois "qu'un électron s'est approprié un proton", c-à-d que sa position tend à se confondre avec celle d'un proton-racine, la couple proton/électron disparaît des radars pour les autres électrons (le champ résultant est nul), et c'est comme si la factorisation se faisait toute seule, sans perte de précision !

Ouais c'est génial. Et aussi, le nombre d'électrons qui vont vers une racine n'est-il pas égal à sa multiplicité ? 

Boécien

Saisissant! Comment se fait-il que ce soit si peu connu?

i.zitoussi

Calli, je ne peux pas te répondre: bien que j'ai mis les cases à cocher dans l'interface, je ne l'ai pas implémenté derrière. Je vais faire des essais. Mais la page wiki et l'article indiquent que oui, et que la convergence devient linéaire au lieu de cubique, donc moins bonne.

Boécien : sans en être certain, je dirais que c'est probablement très connu dans les milieux concernés. Mais c'est vrai que ça mérite d’être connu aussi ailleurs, tellement c'est sexy.

Rescassol

Bonjour
J'ai voulu faire quelques essais, ce n'est pas très concluant.
Voilà mon code en Python (il est possible qu'il y ait des erreurs) et une figure d'un cas où ça converge.
J'ai choisi l'exemple de polynôme de l'article de Aberth et son initialisation avec un peu de hasard en plus.
Les commentaires contiennent des tentatives d'initialisation tirées de l'article de Bini cité dans Wikipédia.

###################################################################
# Résolution d'équations polynômiales par l'algorithme d'Aberth
# Rescassol - 05 Février 2022
###################################################################

###################################################################
# Importations
###################################################################

import numpy as np
from numpy.random import rand
import matplotlib.pyplot as plt

###################################################################
# Fonctions
###################################################################

def maxfig():
    mng = plt.get_current_fig_manager()
    mng.window.showMaximized()

def Tracer_Cercle(xc,yc,R,couleur='r',epaisseur=2):
    t=np.arange(0.0,2*np.pi,0.01)
    x=xc+R*np.cos(t)
    y=yc+R*np.sin(t)
    plt.plot(x,y,couleur,linewidth=epaisseur)

###################################################################
# Script principal de test
###################################################################

# Exemple de polynôme P(z)= z^5 - 10z^4 + 43z^3 - 104z^2 + 150z -100
P=[1,-10,43,-104,150,-100] # Racines 1 +/- 2i , 2 ; 3 +/- i
Deg=5 # Degré de P
D = np.arange(Deg) # Entiers 0 ... Deg-1

# U0=1/(1+np.max([abs(P[i]/P[0]) for i in D+1]))
# Vn=1+np.max([abs(P[i]/P[Deg]) for i in range(Deg)])
# U=[np.max([abs(P[i]/P[k])**(1/(k-i)) for i in range(k)]) for k in D+1]
# V=[np.min([abs(P[i]/P[k])**(1/(k-i)) for i in range(k+1,Deg+1)]) for k in D]
# U=np.concatenate(([U0],U))
# V=np.concatenate((V,[Vn]))

###################################################################

plt.close()
plt.figure(figsize = (20,10), facecolor = 'gold',edgecolor = 'darkgreen')
plt.title('Algorithme de Aberth',size=40,color='m',fontweight='bold')

R, a, C = 3.0, np.pi/(2*Deg), -P[1]/Deg
Z = C + R*np.exp(a*(4*D+1)*1j+np.pi*rand()*1j)#*rand(Deg)
# Z = U[1:]*np.exp(2*np.pi*rand(Deg)*1j)#*rand(Deg)

# print('Z= ',Z.round(4))
Tracer_Cercle(C.real,C.imag,R)
plt.plot(Z.real,Z.imag,'or',markersize=5)

Nmax=10
for n in range(Nmax):
    S=np.zeros(Deg,'complex')
    for k in D :
        Zk=np.concatenate((Z[:k],Z[k+1:]))
        S[k]=sum(1/(Zk-Z[k]))
    T=np.polyval(P,Z)/np.polyval(np.polyder(P),Z)
    W=T/(1-S*T)
    #print('W =',W.round(4))
    for k in D :
        ZZ=np.array([Z[k],Z[k]-W[k]])
        plt.plot(ZZ.real,ZZ.imag,'b')
    Z=Z-W
plt.plot(Z.real,Z.imag,'om',markersize=5)
Rac=np.roots(P)
plt.plot(Rac.real,Rac.imag,'sr',markersize=10)

print('Z =',Z.round(4))
Erreur=np.abs(np.polyval(P,Z))
print('\nErreur =',Erreur)

plt.axis('equal')
plt.show()

Réponse :

Z = [3.+1.j 1.+2.j 2.-0.j 1.-2.j 3.-1.j]

Erreur = [2.92964275e-14 6.55088034e-14 8.06196122e-06 6.35528743e-14 1.76673847e-13]

Cordialement,

Rescassol

i.zitoussi

Bonjour Rescassol,

j'ai aussi fait un petit programme (en html / javascript), et l'algorithme s'avère infaillible.

Ce que je n'ai pas dit, mais que je vais dire tout de suite, c'est que les racines sont connues dès le départ, c'est-à-dire que l'expression $\frac{P'(z)}{P(z)}$ est calculée non pas en calculant séparément $P(z)$ et $P'(z)$, mais directement par $\sum_i\frac{1}{z-z_i}$, où les $z_i$ sont les racines, connues.

C'est d'une part pour éviter le problème des très grands nombres (pas certain que javascript soit fait pour ça, je n'y connais rien), mais aussi parce que je voulais voir le comportement de l'algorithme pour des configurations de points données à l'avance.

Il s'avère que ça converge toujours.

Et pour répondre en même temps à Calli:

Dans la première copie d'écran, il y 6 protonoïdes de charge +100 chacun disposés en arc de cercle, en rouge, (donc le polynôme a 6 racines distinctes de multiplicité 100 chacune), et 300 électrons (points-tests) disposés sur deux cercles, en vert: les protons peuvent accueillir plusieurs électrons, mais la convergence est extrêmement lente.

Dans la deuxième, très "dézoomée", il y a 1 protonoïde (charge +100) et 117 électrons. Pareil, la convergence est désespèrement lente. Il y a du grabuge au début, et 17 électrons ont le bon goût de se désister. Ou se font éjecter. Pas évident de coprendre la psychologie des électrons.

i.zitoussi

Pour la Saint-Valentin.

Calli

Merci pour ton complément de réponse @i.zitoussi.