Wikipedia : Buddhabrot (FR), Buddhabrot (EN)
Flickr Hive Mind : The world's best photos of Buddhabrot
Cette page poursuit l'étude des fractales obtenues par itération d'une fonction complexe, et le cas particulier des ensembles de Mandelbrot et Julia engendrés par les fonctions quadratiques, en s'intéressant à un algorithme particulier de représentation.
Comment se comportent les orbites à l'intérieur de l'ensemble de Mandelbrot, que l'on colorie souvent en noir uniforme ? On sait que les orbites restent à l'intérieur de l'ensemble (par définition), mais certains points sont plus souvent visités que d'autres.
Pour mieux observer et comprendre cela, on peut considérer un grand nombre d'orbites (pour un grand nombre de points de départ), et colorier chaque point du plan en fonction du nombre de fois où il a été visité par des orbites. On obtient ainsi une densité de passage des orbites au voisinage de chaque point.
De nombreuses variantes existent ; la plus célèbre consiste à tracer les orbites des points qui ne sont pas dans l'ensemble de Mandelbrot. Dans ce cas, l'image obtenue est appelée Buddhabrot (voir The Buddhabrot technique par Melinda Green), mais remarquons qu'elle nécessite typiquement des centaines de milliers d'itérations pour donner un résultat intéressant. La méthode choisie ici (densité des orbites des points appartenant à l'ensemble de Mandelbrot) est "efficace" avec un petit nombre d'itérations, et l'image obtenue est parfois appelée anti-Buddhabrot.
Ci-dessous, le code Python permettant de générer cette image (en plusieurs centaines de secondes).
""" Fractal orbit density map plotter ('Buddhabrot') creates a png image """ import datetime timedebut = datetime.datetime.now() print"\n Fractal orbit density map plotter \n" # from numpy import sqrt, empty, real, imag, pi, cos, sin, exp, sign, array, max, log, zeros, random from numpy import empty, real, imag, pi, cos, sin, log, zeros, max from cmath import polar, rect from PIL import Image # Some parameters # Max number of iterations maxiter = 100 # resolution finale 2D imgx = 576 # imgy = 1016 imgy = 508 # size of the cloud of points ptx = 1200 pty = 800 numpts = ptx * pty # initialisation cr = empty(imgx * imgy) cr.shape = imgx, imgy cg = empty(imgx * imgy) cg.shape = imgx, imgy cb = empty(imgx * imgy) cb.shape = imgx, imgy cr=zeros(imgx * imgy) cg=zeros(imgx * imgy) cb=zeros(imgx * imgy) cr.shape = imgx, imgy cb.shape = imgx, imgy cg.shape = imgx, imgy z = empty(maxiter + 1, dtype = 'complex') nuage = empty(numpts, dtype = 'complex') image = Image.new("RGB", (imgy, imgx)) # initial nuage construction pt = 0 for i in range(ptx): for j in range(pty) : nuage[pt] = -2.5 + i / float(ptx) * 5.0 + 1j *( -2.5 + j / float(pty) * 5.0 ) pt += 1 # for each point in the initial nuage... for pt in range(numpts) : z1 = nuage[pt] # color depending on the starting point x0 = (real(z1) + 2.0) / 3.0 y0 = (imag(z1) + 1.25) / 2.5 tb = .5 * ( cos( pi * x0) + 1.0) tr = 1.0 - tb tg = .5 * ( cos( pi * y0) + 1.0) # fractal parameters (Mandelbrot/Julia) c = z1 z[0] = z1 # ... we compute the orbit under iteration of the quadratic polynomial iter = 0 while ( iter < maxiter ) and ( abs(z[iter]) < 5 ) : iter += 1 z[iter] = z[iter-1] * z[iter-1] + c # and if the starting point belongs to the Mandelbrot, we count (plot) each point of the orbit if iter == maxiter : for iter in range(maxiter / 2, iter + 1) : w = z[iter] x = real(w) y = imag(w) kx = imgx*(x + 2.0) / 3.0 ky = imgy*(y + 1.25) / 2.5 if (kx >= 0) and (kx < imgx) and (ky >= 0) and (ky < imgy) : cg[kx,ky] += tg cb[kx,ky] += tb cr[kx,ky] += tr # normalization coefficients limr = log(max(cr) + 1.0) if limr == 0.0 : limr = 1.0 limg = log(max(cg) + 1.0) if limg == 0.0 : limg = 1.0 limb = log(max(cb) + 1.0) if limb == 0.0 : limb = 1.0 # plot the normalized logarithm of the density (inversion x/y) for jx in range(imgx): for jy in range(imgy): coulr = 255 * log(cr[jx, jy] + 1) / (limr) coulg = 255 * log(cg[jx, jy] + 1) / (limg) coulb = 255 * log(cb[jx, jy] + 1) / (limb) image.putpixel((jy , jx), (coulr, coulg, coulb)) image.save('buddhabrot.png', "PNG") timefin = datetime.datetime.now() print print "Total execution time : ", (timefin-timedebut).seconds, "s"
Pour démarrer en Python, installez la distribution Python(x,y).