Densité des orbites sous l'itération d'une fonction - Le "Buddhabrot"

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.

buddhabrot.png

Ci-dessous, le code Python permettant de générer cette image (en plusieurs centaines de secondes).

python-logo.gif
"""
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).

Unless otherwise stated, the content of this page is licensed under Creative Commons Attribution-ShareAlike 3.0 License