MMI, Lyon
Il est temps de faire un peu de maths.
On va considérer la fonction f(z)=exp(z−3/2) où z est une variable complexe.
Elle définit un système dynamique, une équation d'évolution à temps discret. L'état au temps n est donné par un nombre complexe zn et l'évolution par
L'ensemble de Julia, noté J(f) ou simplement J, a plusieurs définitions qui donnent le même objet. Je ne vais pas les rappeler ici. En voici une, volontairement vague : c'est le lieu où la dynamique est instable. Je donne ci-dessous un moyen pratique de déterminer J dans notre cas particulier.
La fonction f choisie ci-dessus possède un point fixe attractif a=0.3107... Le bassin d'attraction est défini comme l'ensemble des points dont l'orbite est attirée par a. Ici, le bassin est connexe, non borné, et simplement connexe (pas de trou), et enfin : l'ensemble de Julia est le complémentaire du bassin.
Cela suggère un premier algorithme pour tenter de dessiner J(f). Chaque pixel de l'image correspond à une valeur de z0. On demande à l'ordinateur de calculer la suite zn. Si un des zn tombe assez près de a, on colorie de pixel en blanc. Si n atteint une valeur trop élevée, on le colorie en noir. J'ai écrit tenter car on ne sais pas à ce stade si ça va bien fonctionner. Par exemple il y a un risque que J se faufile entre les valeurs z0 que l'on échantillonne et qu'on obtienne essentiellement du blanc.
Voici le programme.
Au lieu d'un test de proximité au point fixe attractif 'a', on effectue le test "Re(z)<1". Il se trouve en effet que le demi-plan "Re(z)<1" est inclus dans le bassin d'attraction et contient 'a'. Son image est la boule de centre 0 et de rayon 0.606..., épointée.
from __future__ import division import png from array import array import cmath center = 4 + 0j xw = 12.0 width = 600 height = 600 max_iter = 50 NONE=0 ; MAX=1 ; BASIN=2 def compute(j,i) : z = center + (j-width/2+(height/2-i)*1.0j)*xw/width k = 0 reason = NONE while True: if k>=max_iter : reason = MAX; break if z.real < 1. : reason = BASIN; break z = cmath.exp(z-1.5) k += 1 if reason==BASIN : return [255,255,255] return [0,0,0] pixels=[] for i in xrange(height) : row = [] for j in xrange(width) : row.extend(compute(j,i)) pixels.append(row) image = png.from_array(pixels,'RGB') image.save('essai-2.1.png')En vert, les parties nouvelles/modifiées par rapport à la version précédente du programme. Les paramètres center, xw, width, height contrôlent la géométrie de la fenêtre. max_iter règle le nombre maximal d'itérations de f pour un pixel donné.
from __future__ import division permet la chose suivante : que la division de deux variables de type integer donne une variable de type float. Traditionnellement, ça retournerait un entier : le quotient de la division euclidienne, ce qui n'est pas très pratique (dans 90% des cas j'ai besoin du quotient réel).
Python intègre les nombres complexes, de base ! Un nombre complexe s'écrit par exemple 0.2-0.3j avec j2 = −1. Par contre il faut faire import cmath pour avoir accès à certaines variantes complexes de fonctions du type exp, cos, etc...
Une boucle while condition : exécute le bloc qui suit en boucle tant que la condition donne "True". La commande break permet de sortir immédiatement de la boucle la plus interne contenant le break (même si celui est dans un sous-bloc d'un autre type, "if" dans notre exemple). Le pointeur d'exécution se retrouve alors juste après le bloc de cette boucle.
Note. J'ai lu il y a longtemps qu'il fallait éviter les "break" et encore pire les "goto". L'argument est que cela complique le déboggage, ce qui est vrai. Dans mon cas cependant, ces pratiques rendent le code à la fois plus lisible et plus facile à écrire. J'ai essayé pendant plusieurs années les 2 méthodes, sans et avec, j'ai donc bien pu comparer. Mes programmes en C++ qui calculent des images sont courts et incluent souvent des "goto". Python n'a pas de commande "goto" mais a le "break".
On lance le programme depuis un terminal et on obtient... on n'obtient pas d'image mais un message :
Traceback (most recent call last):
File "prog-2a.py", line 31, in <module>
row.extend(compute(j,i))
File "prog-2a.py", line 21, in compute
z = cmath.exp(z-1.5)
OverflowError: math range error
C'est un message d'erreur, avec trop d'information ce qui le rend difficile à lire. L'important est à la fin : c'est une erreur d'overflow, c.à.d. on a une opération qui retourne un réel trop grand par rapport à ce que permet la représentation machine : les flottants. Il semblerait que cette représentation dépende de la machine exécutant le programme Python. En général aujourd'hui un flottant est représenté sur 64 bits : un bit de signe, un bit pour le signe de l'exposant, 10 bits pour l'exposant, le reste pour la mantisse, soit 52 bits. Ainsi le flottant maximal est proche de 21024 ≈ 1.8E308.
Comme on l'a vu, en cas d'overflow, une erreur a lieu et par défaut Python s'interrompt. On aimerait pouvoir continuer le flot, par exemple en mettant le pixel en rouge et en passant au suivant. C'est possible: nous verrons plus tard un moyen d'intercepter l'erreur et d'éviter l'interruption. Mais pour l'instant, on va essayer d'éviter les dépassements.
Comme le module de exp(z-1.5) vaut exp(Re(z)-1.5), l'overflow a lieu quand on applique f à un nombre z dont la partie réelle est trop grande. On introduit donc ci-dessous un nouveau paramètre seuil, partie réelle à droite de laquelle on va stopper l'itération. Pour atteindre le flottant maximal avec exp(x), il faut que x vaille environ 700, nous allons donc prendre un seuil inférieur à cette valeur.
Quelle couleur donner au pixel quand l'orbite dépasse le seuil ? On considère qu'un point qui se retrouve suffisamment à droite est forcément parti proche de J(f), ce qui dans notre cas est une assez bonne approximation. (On justifiera ça plus loin et on donnera également des améliorations).
Ci-dessous, le programme intégrant le test de seuil, les ajouts sont indiqués en rouge. J'ai changé le choix de couleurs des pixels : c'est toujours en blanc pour ceux qui passent dans le demi-plan Re(z)<1 avant max_iter mais il y a maintenant des points bleus : ceux qui finissent par vérifier Re(z)>seuil en moins de max_iter itérations, enfin je colorie en rouge (au lieu de noir) ceux qui ont atteint max_iter avant de vérifier l'une des deux conditions.
from __future__ import division import png from array import array import cmath center = 4 + 0j xw = 12.0 width = 600 height = 600 max_iter = 50 seuil = 40.0 # more than 700 -> overflow error NONE=0 ; MAX=1 ; BASIN=2 ; SEUIL=3 def compute(j,i) : z = center + (j-width/2+(height/2-i)*1.0j)*xw/width k = 0 reason = NONE while True: if k>=max_iter : reason = MAX; break if z.real > seuil : reason = SEUIL; break if z.real < 1. : reason = BASIN; break z = cmath.exp(z-1.5) k += 1 if reason==BASIN : return [255,255,255] if reason==SEUIL : return [0,0,255] return [255,0,0] pixels=[] for i in xrange(height) : row = [] for j in xrange(width) : row.extend(compute(j,i)) pixels.append(row) image = png.from_array(pixels,'RGB') image.save('essai-2.2.png')Voici le résultat avec seuil=40.
Pour comparer, j'ai fait tourner le programme avec seuil=500.
On observe des moirés conséquents. Il est intéressant de noter l'absence de points rouge : aucun des 6002=360 000 points de l'échantillon n'a atteint 50 itération avant que la partie réelle de l'orbite ne sorte de l'intervalle [1,500].