TangenteX.com

Dessiner l'ensemble de Mandelbrot

Nous allons aborder dans cette page deux outils fondamentaux utilisés en physique : les nombres complexes et plus particulièrement les suites de nombres complexes et les fractales, objets mathématiques exotiques dont les dimensions ne sont pas entières.

Tous mes lecteurs ont sans aucun doute déjà rencontré les nombres complexes. Nous les utilisons pour travailler sur des fonctions périodiques parce qu'ils facilitent grandement les calculs. Nous les utilisons aussi en physique quantique, ses équations sont toutes des équations dans le corps des complexes \( \mathbb{C}\). Mais il est moins courant d'utiliser des suites complexes, comme nous allons le faire ici.

Les fractales sont des objets mathématiques qui servent à modéliser bon nombre de systèmes physiques. Par exemple, on les rencontre quand on étudie les problèmes de diffusion. Les informaticiens les connaissent surtout pour leur capacité à créer des représentations réalistes d'objets comme des arbres, des réseaux biologiques ou encore des reliefs géographiques. Wikipedia vous en dira beaucoup plus ici.

Dans cette page estivale, m'est venue l'envie de construire une figure très célèbre chez les mathématiciens, physiciens et informaticiens. Il s'agit de l'ensemble de Mandelbrot. Vous trouverez sur le net une multitude de représentations de cet ensemble, mais avouez qu'il est plus amusant de le construire soi-même, non ?

Un peu de mathématique

L'ensemble de Mandelbrot

Ce que nous appelons 'ensemble de Mandelbrot' n'a pas été inventé (ou découvert, comme vous voudrez) et étudié par Benoît Mandelbrot (X44) mais par deux mathématiciens français du début du XXème siècle : Gaston Julia et Pierre Fatou. Mandelbrot fut le premier à simuler sur ordinateur cet ensemble au début des années 80 et à en créer des représentations graphiques. Enfin, la définition moderne de cet ensemble et son nom sont dûs à un autre mathématicien français, Adrien Douady.

L'ensemble de Mandelbrot est défini à partir de la relation de récurrence sur  \( \mathbb{C}\) :

\( \left\{ \begin{array}{lr} z_{n+1} = z_n^2 + c \\ z_0 = 0  \\ \end{array} \right. \)

avec c une constante complexe quelconque et z une variable complexe. Pour information, on appelle 'orbite' la suite des nombres complexes obtenus à partie de c.

L'ensemble de Mandelbrot \( \mathbb{M}\) est l'ensemble des nombres c de \( \mathbb{C}\) tel que la suite \( {z_n} \) soit bornée. Si la suite \( {z_n} \) est bornée en c, alors c est un élement de \( \mathbb{M} \), sinon c n'appartient pas à \( \mathbb{M} \). C'est aussi simple que ça...

Rappelons que la suite \( {z_n} \) est bornée si la valeur du module de \( z_n \), quand n tend vers l'infini, est toujours plus petite ou égale qu'un nombre réel fixé, nombre que l'on nomme traditionnellement l'horizon de la suite. Si la suite n'est pas bornée, on dit qu'elle diverge, c'est à dire que lorsque n tend vers l'infini, le module de \( z_n \) tend vers l'infini.

Comment le construire

Commençons par délimiter le domaine du plan complexe sur lequel nous allons travailler. Puis pour chaque point de ce domaine, calculons la suite \( {z_n} \) pour des valeurs de n comprises en 0 et un nombre "raisonnable" d'itérations. La valeur de "raisonnable" dépend de nos moyens ! En théorie, il faudrait faire tendre n vers l'infini, mais comme l'infini n'existe pas en informatique, on dira que n = 100 est raisonnable. C'est arbitraire, et vous pourrez faire des essais avec une autre valeur. Sachez que c'est surtout la puissance de votre processeur qui vous limitera...

En chaque point de notre domaine, nous allons déterminer si la suite est bornée, en comparant le module de \( z_n \) à l'horizon. Si la suite est bornée, le point appartient à \( \mathbb{M}\) et nous le dessinerons en noir. Sinon, la couleur du point déterminera la rapidité de la divergence de la suite.

Une première approche

L'algorithme étant particulièrement simple, je vais directement passer à la programmation. Celle des fonctions de calcul ne pose aucun problème. Et matplotlib offre les outils nécessaires pour visualiser de façon satisfaisante les résultats.

Dans une première approche, je ne vais pas m'occuper des performances du code, mais coder simplement l'algorithme décit au chapitre précédent.

Définir les paramètres de calcul et d'affichage

Le domaine exploré

Les bornes du domaine que j'ai choisi d'explorer sont définies dans le code :

xmin, xmax = -2.25, +0.5
ymin, ymax = -1.25, +1.25

Les limites sont choisies en faisant un rapide calcul sur le domaine intéressant à explorer pour obtenir une vue globale de \( \mathbb{M}\). Elles peuvent être adaptées, en particulier si vous voulez zoomer sur une zone particulière du domaine. Nous verrons cela plus loin.

Les paramètres de calcul

Les voici :

nbiter = 100
horizon = 100.0

avec :

Les paramètres d'affichage

J'ai choisi les valeurs de définition de l'affichage suivantes :

LX = 2000
LY = 1500
width =  15
height = 15*LY/LX

Le choix de ces valeurs dépend de votre écran (je travaille sur un Mac avec un écran Rétina). Attention, elles influent grandement sur le temps de calcul (elles fixent le pas des boucles que nous allons voir ci-dessous) !

Les fonctions de calcul

Calcul de la suite zn pour un point

Le code de cette fonction est conforme à l'algorithme décrit ci-dessus. Les commentaires dans le programme sont, je crois, suffisants à sa compréhension.

def MandelbrotDef(c,niter,horizon): 
    z = c                     # initialisation de la récurrence
    for n in range(niter):    # déroulement du nombre d'itérations max
        if abs(z) > horizon:  # test de divergence
            return n          # si divergence, on retourne le nb d'itérations avant divergence     
        z = z**2 + c          # calcul du terme suivant de la suite
    return 0                  # la suite est convergente, on retourne 0

Construction de l'ensemble de Mandelbrot

Je définis deux vecteurs qui contiennent les coordonnées x et y, soit la partie réelle et imaginaire de chaque affixe du domaine complexe. Les caractéristiques de ces vecteurs sont passées en paramètres de la fonction. Ils ont été définis ci-dessus. Je construis aussi le vecteur P réel qui me servira à stocker le nombre d'itérations effectués à chaque point.

Puis je parcours le domaine complexe avec deux boucles imbriquées, en lançant à chaque point le calcul de la suite. Ce qui me donne le code suivant :

def MandelbrotSet(xmin, xmax, ymin, ymax, Lx, Ly, niter):
    r1 = linspace(xmin, xmax, Lx)
    r2 = linspace(ymin, ymax, Ly)
    P = zeros((Lx,Ly))
    for i in range(Lx):
        for j in range(Ly):
            P[i,j] = MandelbrotDef(r1[i] + 1j*r2[j], niter)
    return P

Le calcul

Rien à signaler, si ce n'est la mesure du temps de calcul et son affichage à la fin du calcul :

debut = time.time()
Points = MandelbrotSet(xmin, xmax, ymin, ymax, LX, LY, nbiter, horizon)
print("Temps de calcul : %s s" % (time.time() - debut))

Le code d'affichage

J'utilise la fonction imshow() de matplotlib, bien adaptée à cet usage :

# préparation de l'affichage
fig = plt.figure(figsize=(width, height))
fig.suptitle("Ensemble de Mandelbrot", fontsize=16, color = 'w')
ax = fig.add_axes([0.0, 0.0, 1.0, 1.0], frameon=False, aspect=1)
# paramètres de réglage de l'affichage et affichage
gamma = 0.4
cmap = 'hot'
plt.imshow(Points.T, extent=[xmin, xmax, ymin, ymax], interpolation="bicubic",
           cmap = cmap, norm = colors.PowerNorm(gamma))
plt.show()

Une petite explication à propos de la fonction PowerNorm(). Elle permet de régler l'équilibre et l'intensité des couleurs. Je trouvais que la figure brute était trop sombre. Pour régler l'intensité et l'équilibre des couleurs, modifiez la valeur du paramètre gamma.

Si vous le souhaitez, vous pouvez utiliser un autre jeu de couleurs (paramètre cmap). J'ai utilisé le jeu 'hot' par classicisme, mais vous pouvez tout à fait en choisir un autre.

Le code complet du script est disponible ici.

Le résultat

Le temps de calcul est relativement long... Avec 100 itérations par point sur le domaine et pour les paramètres d'affichage choisis, j'obtiens un temps de calcul de 96 s environ.

Voici la figure obtenue :

Tous les points du domaine complexe qui sont en noirs sont des éléments de \( \mathbb{M}\). Pour les autres points, la couleur varie du jaune très claire, pour les points où la suite diverge lentement, au rouge sombre, pour les points où la suite diverge rapidement. Plus un point est éloigné du centre de \( \mathbb{M}\), plus la divergence de la suite en ce point est rapide.

Sur cette figure, vous pouvez distinguer différentes formes :

Ces formes ont donné lieu à de nombreuses études mathématiques. Vous pouvez calculer assez facilement le centre de la cardioïde et des bourgeons ainsi que les extrémités des antennes. Il suffit (!) de traiter la suite complexe comme un polynôme complexe et de chercher ses racines.

Un peu d'optimisation

Vous avez noté que le temps de calcul du script précédent est assez long. La faute à la double boucle dans la fonction MandelbrotSet() et aussi à la méthode de détermination de la divergence. D'aucuns vous diront qu'il faut virer le calcul du module (la fonction abs()) qui utilise une racine carrée et faire le test sur le carré du module, mais c'est assez peu efficace ! La fonction abs() est une fonction standard de Python et elle est déjà optimisée...

Il est bien préférable d'utiliser les fonctions natives de numpy (ou de scipy, au choix), qui sont codées en C et déjà optimisées. Et aussi de vectoriser le plus possible...

Je vous propose donc le code suivant dans la fonction MandelbrotSetSpeed(), qui remplace les fonctions MandelbrotDef() et MandelbrotSet().

def MandelbrotSetSpeed(xmin, xmax, ymin, ymax, Lx, Ly, niter, horizon):    
    r1 = linspace(xmin, xmax, Lx)
    r2 = linspace(ymin, ymax, Ly)
    C = r1 + r2[:, None]*1j
    P = zeros(C.shape, dtype=int)
    Z = zeros(C.shape, complex64)
    for n in range(niter):
        I = less(abs(Z), horizon)
        P[I] = n
        Z[I] = Z[I]**2 + C[I]
    P[P == niter-1] = 0
    return P

On retrouve dans ce code la même structure que dans le précédent. Seules diffèrent l'utilisation de la vectorisation (pour c) et l'usage de la fonction less() pour comparer le module de z et l'horizon.

Pour le reste, le code du script est identique au script précédent, il est disponible ici.

Le résultat

Le temps de calcul est beaucoup plus court... Avec 100 itérations par point sur le domaine et pour les paramètres d'affichage choisis, j'obtiens un temps de calcul de 5 s environ, soit un code presque 20 fois plus rapide !

Voici la figure obtenue, identique à la précédente :


Il existe d'autres moyens plus ou moins efficaces pour optimiser ce code en Python. Si vraiment vous voulez un code rapide, je vous conseille d'écrire votre programme en C et d'utiliser gnuplot pour l'affichage. Avec les paramètres du script, vous devriez obtenir un temps de calcul inférieur à la seconde.

Des pistes d'approfondissement

A propos de l'ensemble de Mandelbrot

Wikipedia propose une étude détaillée de l'ensemble de Mandelbrot. Vous trouverez aussi sur le net des analyses mathématiques plus poussées sur les fractales, les ensembles de Julia et leurs applications en physique. Voir aussi ici une approche sympathique de l'ensemble de Mandelbrot.

Une autre approche intéressante est de modifier la relation de récurrence qui définit l'ensemble. Essayez par exemple la suite :

\( \left\{ \begin{array}{lr} z_{n+1} = z_n^3 + c \\ z_0 = 0  \\ \end{array} \right. \)

ou encore la suite :

\( \left\{ \begin{array}{lr} z_{n+1} = e^z + c \\ z_0 = 0  \\ \end{array} \right. \)

ou toute autre suite complexe que vous pourriez imaginer...

Vous pouvez aussi, ce qui vous sera profitable pour vos études, étudier de manière approfondie les suites complexes et les fonctions holomorphes (les fonctions complexes). C'est un sujet d'étude et de réflexion insondable. Il faut avouer que les physiciens voient plutôt ça comme des outils que comme des sujets de réflexion, mais bon..

Sur le dessin de l'ensemble de Mandelbrot

Faire des zooms

Avec les paramètres suivants :

xmin, xmax = -0.8, -0.7
ymin, ymax = 0.0, 0.1
nbiter = 1024

voilà l'image que j'obtiens :


Autre essai, avec les paramètres suivants :

xmin, xmax = -0.755, -0.740
ymin, ymax = 0.06, 0.07
nbiter = 2048

ce qui me donne :

Le temps de calcul reste raisonnable : 120 secondes environ sur mon Mac I5.

En piochant dans les détails d'une image lors d'un essai, avec les paramètres suivants:

xmin, xmax = -0.748770, -0.748730
ymin, ymax = 0.06505, 0.06510
nbiter = 2048

j'ai obtenu l'image suivante :

l'image d'un ensemble de Mandelbrot ! Nous sommes bien dans une structure fractale !

Lisser les couleurs représentant la rapidité de divergence

Dans les figures globales de l'ensemble de Mandelbrot, vous avez constaté que le dégradé des couleurs à l'extérieur de l'ensemble n'était pas terrible. Il est possible d'obtenir une transition plus douce, et même un lissage (anti-aliasing) de ces couleurs.

Le principe est proposé par Linas Vepstas. Il s'agit d'utiliser un horizon beaucoup plus, 240 selon l'auteur, et procéder à un calcul particulier sur les couleurs. La justitication théorique du procédé est fournie ici, par Linas Vepstas.

J'ai modifié le code de mon algorihme naïf pour ce traitement :

def MandelbrotDef(c,niter, horizon, logh): 
    z = c                    
    for n in range(niter):   
        if abs(z) > horizon: 
           return n - log(log(abs(z)))/log(2) + logh     
        z = z**2 + c          
    return 0                  

et la définition de l'horizon :

horizon = 2.0**40
logh = log(log(horizon))/log(2)

J'obtiens la figure suivante :

Vous constatez que le dégradé des couleurs est bien plus agréable. Mais le temps de calcul est aussi bien plus long, environ 260 secondes sur mon Mac.

Le code complet du script est ici.

Contenu et design par Dominique Lefebvre - www.tangenteX.com août 2019 Licence Creative Commons Contact :