Mandelbrotova množica

preizkusi program

V tem članku bom čim bolj preprosto razložil mandelbrotovo množico iz matematičnega vidika in pokazal, kako se naredi preprost program za vizualizacijo mandelbrotove množice.

Mandelbrotova množica je množica kompleksnih števil \(c\), pri katerih funkcija \(f(z)=z^2+c\) ne narašča v neskončnost (vir: https://en.wikipedia.org/wiki/Mandelbrot_set). To funkcijo lahko zapišemo tudi kot zaporedje \(z_{i}=z_{i-1}^2+c\). Pri tem je \(z\) kompleksno število, ki je na začetku \(z=0+0i\), potem pa se pri vsaki ponovitvi spremeni. Velika sprememba pomeni, da \(z\) narašča v neskončnost, in da \(c\) ni del mandelbrotove množice.

Množico kompleksnih števil \(c\) lahko prikažemo v koordinatnem sistemu tako, da realni del prikažemo na \(x\), imaginarni del pa na \(y\) osi. Ko pobarvamo vse, ki spadajo v mandelbrotovo množico, vidimo zanimive fraktale.

Poglejmo si primer

Na sliki je prikazano število \(c=2+3i\).

Ugotovimo, če to število spada v mandelbrotovo množico:

V prvem koraku izračunamo \(z_1\), tako da ga nastavimo na vrednost izraza \(z_0^2+c\) kjer je \(z_0=0+0i\) in \(c=2+3i\):

\(z_1=z_0^2+c\) \(=(0+0i)^2+2+3i=2+3i\)

V drugem koraku, na podoben način izračunamo \(z_2\):

\(z_2=z_1^2+c\) \(=(2+3i)^2+2+3i\) \(=2^2+2*2*3i*(3i)^2+2+3i\) \(=4+12i-9+2+3i\) \(=-3+15i\)

Po enakem postopku izračunamo še nekaj \(z\)-jev (od zdaj naprej bom izpuščal vmesne korake, ker postajajo predolgi):

\(z_3=z_2^2+c\) \(=(-3+15i)^2+2+3i\) \(=-124-87i\)

\(z_4=z_3^2+c\) \(=(-124+87i)^2+2+3i\) \(=38229+37239i\)

\(z_5=z_4^2+c\) \(=(38229+37239i)^2+2+3i\) \(=74713322+2847219465i\)

Vsakič ko izračunamo nov \(z\), naredimo eno ponovitev ali iteracijo.

Že pri petih iteracijah ugotovimo, da se \(z\) zelo hitro povečuje, zato lahko sklepamo, da \(z\) narašča v neskončnost in da \(c=2+3i\) ni v mandelbrotovi množici. Teoretično bi morali izračunati neskončno iteracij, da bi se prepričali, da \(c\) res ni v mandelbrotovi množici, ampak za naš približek je dovolj, da naredimo samo 5 ponovitev. Ko bomo hoteli mandelbrotovo množico izrisati na zaslon, bomo število iteracij povečali.

Poglejmo še kako izgleda iteracija za število, ki je v mandelbrotovi množici. Vzemimo \(c= -1+0.2i\):

\(z_1=z_0^2+c\) \(=(0+0i)^2-1+0.2i\) \(=-1+0.2i\)

\(z_2=z_1^2+c\) \(=(-1+0.2i)^2-1+0.2i\) \(=-0.04-0.2i\)

\(z_3=z_2^2+c\) \(=(-0.04-0.2i)^2-1+0.2i\) \(=-1.0384+0.216i\)

\(z_4=z_3^2+c\) \(=(-1.0384+0.216i)^2-1+0.2i\) \(=0.0316+-0.2485i\)

\(z_5=z_4^2+c\) \(=(0.0316+-0.2485i)^2-1+0.2i\) \(=-1.0607+0.1842i\)

Ko naredimo nekaj iteracij vidimo, da \(z\) ne narašča v neskončnost, zato lahko sklepamo, da \(c\) je v mandelbrotovi množici.

Začnimo programirati

Naredimo program, ki bo ugotovil, ali je število \(c\) del mandelbrotove množice. Za programiranje bom uporabil programski jezik JavaScript.

function mandelbrot(ca, cb){

}

Začnimo z deklaracijo funkcije mandelbrot(), ki za parameter vzame kompleksno število \(c\). Ker pa računalnik lahko računa samo z realnimi števili, sta to v resnici dva parametra ca in cb, ki predstavljata realni in imaginarni del.

var maxIteracij = 50;
var prag = 100;
function mandelbrot(ca, cb){
    var za = 0;
    var zb = 0;
    var n = 0;
}

V funkciji definiramo še kompleksno število \(z\), ki je spet razdeljeno na realni in imaginarni del za in zb in spremenljivko n, ki se bo vsako iteracijo povečala. Izven funkcije pa definiramo dve globalni spremenljivki: maxIteracij, ki določa pri kateri iteraciji se program ustavi, in prag, ki določa kdaj je število dovolj veliko, da lahko rečemo, da narašča proti neskončnosti in se program lahko ustavi.

while(/*je z manjši od prag*/){
    /*z = z^2+c*/
}

V funkcijo dodamo while zanko, ki se ponavlja, dokler je pogoj z<prag resničen. Tukaj spet nastopi težava, ker je \(z\) kompleksno število. Da ga primerjamo z realnim številom prag, lahko uporabimo pitagorov izrek.

Velikost števila \(z\) lahko izračunamo kot razdaljo točke od izhodišča koordinatnega sistema.

\(z=a+bi\)

\(d(z)=\sqrt{a^2+b^2}\)

\(\sqrt{a^2+b^2}<prag\)

Da se izognemo korenu, ki zahteva veliko procesorske moči, cel izraz kvadriramo:

\(a^2+b^2<prag^2\)

Zdaj lahko pogoj vstavimo v program:

while(Math.pow(za,2)+Math.pow(zb,2) < Math.pow(prag,2)){
    /*z = z^2+c*/
}

Vsakič, ko se while zanka ponovi, mora program izračunati eno element zaporedja \(z_{i}=z_{i-1}^2+c\). Ker računalnik računa samo z realnimi števili, moramo izraz s kompleksnimi števili pretvoriti v dva z realnimi.

\(z_{nov}=z^2+c\)

Kompleksni števili \(z\) in \(c\) zapišemo v obliki \(z=a+bi\):

\(z_{nov_a}+z_{nov_b}i\) \(=(z_a+z_bi)^2+c_a+c_b i\)

Razširimo \((z_a+z_b i)^2\) na \(z_a^2+2z_a z_b i-z_b^2\):

\(z_{nov_a}+z_{nov_b}i\) \(=z_a^2+2z_a z_b i-z_b^2+c_a+c_bi\)

Posebej zapišemo vse realne in imaginarne člene (imaginarni so tisti, ki imajo zraven \(i\)):

\(z_{nov_a}=z_a^2-z_b^2+c_a\)

\(z_{nov_b}i=2z_a z_b i+c_b i\)

Pri enačbi za imaginarni del \(z_{nov}\) lahko \(i\) izpostavimo in okrajšamo:

\(z_{nov_b}=2z_a z_b+c_b\)

Zdaj računalnik nima nobenih težav pri reševanju enačbe in jo lahko vpišemo v program:

while(Math.pow(za,2)+Math.pow(zb,2) < Math.pow(prag,2)){
    var zanov = Math.pow(za,2)-Math.pow(zb,2)+ca;
    var zbnov = 2*za*zb + cb;
    za = zanov;
    zb = zbnov;
}

Morda se sprašujete, zakaj potrebujemo spremenljivki zanov in zbnov, zakaj ne bi rezultata nastavili že kar pri izračunu. Če bi za nastavili že pri prvem izračunu, bi pri drugem izračunu računali z napačno vednostjo in bi zb izračunali narobe.

Skoraj smo končali, dodati moramo samo še pogoj, ki bo zagotovil, da tudi pri številih, ki nikoli ne dosežejo praga v razumnem času končamo izvajanje funkcije. To bomo naredili tako, da če števec iteracij n preseže maxIteracij prekinemo izvajanje while zanke.

var maxIteracij = 50;
var prag = 100;
function mandelbrot(ca, cb){
    var za = 0;
    var zb = 0;
    var n = 0;
    while(Math.pow(za,2)+Math.pow(zb,2) < Math.pow(prag,2)){
        var zanov = Math.pow(za,2)-Math.pow(zb,2)+ca;
        var zbnov = 2*za*zb + cb;
        za = zanov;
        zb = zbnov;
        
        n++;
        if(n>=maxIteracij)break;
    }
    return n;
}

Na koncu dodamo ukaz return n, da nam bo funkcija vrnila koliko iteracij je trajalo preden se je \(z\) približal neskončnosti. Če se \(z\) nikoli ni približal neskončnosti, bo funkcija vrnila vrednost maxIteracij.

V funkcijo lahko dodamo še ukaz console.log(), ki bo vsako iteracijo v konzolo izpisal števec iteracij in vrednost števila \(z\).

var maxIteracij = 50;
var prag = 100;
function mandelbrot(ca, cb){
    var za = 0;
    var zb = 0;
    var n = 0;
    while(Math.pow(za,2)+Math.pow(zb,2) < Math.pow(prag,2)){
        var zanov = Math.pow(za,2)-Math.pow(zb,2)+ca;
        var zbnov = 2*za*zb + cb;
        za = zanov;
        zb = zbnov;
        console.log(n+". z="+za+"+"+zb+"i");
        n++;
        if(n>=maxIteracij)break;
    }
    return n;
}

Tega ukaza ne smemo pozabiti odstraniti, ko bomo preverjali več števil, ker zelo upočasni program.

Tukaj lahko preizkusite, program, ki ugotovi ali je poljubno število v mandelbrotovi množici. Ob vsaki ponovitvi izpiše zaporedno število ponovitve in vrednost \(z\)-ja. Če želite, da program naredi več iteracij, povečajte spremenljivki maxIteracij in prag.



            
        

Risanje fraktalov

V nadaljevanju bomo naredili funkcijo, ki bo za vsak piksel na zaslonu preverila, ali je del mandelbrotove množice, in ga ustrezno pobarvala.

Najprej na vrhu programa definirajmo še nekaj spremenljivk:

var c = document.getElementById("canvas");
var ctx = c.getContext("2d"); //objekt s funkcijami za risanje
var w = c.width = 600;	 //širina risalne površine (600)
var h = c.height = 400;    //višina (400)
var zoom = 100; //velikost fraktala
var zamikX = 0; //zamik po x osi
var zamikY = 0; //zamik po y osi

Nato naredimo funkcijo izris(), ki ima dve for zanki, ki gresta čez vse piksle na risalni površini.

function izris(){
    for(var x=0; x<w; x++){
        for(var y=0; y<h; y++){
            
        }
    }
}

V zanki dodamo spremenljivki ca in cb, ki predstavljata kompleksno število \(c\) za katerega izračunamo ali pripada mandelbrotovi množici, tako da kličemo funkcijo mandelbrot(). Glede na rezultat te funkcije pobarvamo piksel (x, y). Spomnimo se, da funkcija mandelbrot() vrne število iteracij pri katerih število \(z\) prekorači prag, če ga ne prekorači nikoli nam funkcija vrne število maxIteracij.

Vsak piksel ima lahko barvo rgb(r,g,b), kjer r pomeni koliko je rdeče, g koliko zelene in b koliko modre barve. Če imajo r, g in b vrednost nič dobimo črno barvo, če pa imajo največjo možno vrednost 255, pa dobimo belo barvo. Da piksel pobarvamo glede pripadnost mandelbrotovi množici, nastavimo r, g in b na vrednost funkcije mandelbrot().

function izris(){
    for(var x=0; x<w; x++){
        for(var y=0; y<h; y++){
            var ca = x;
            var cb = y;
            var v = mandelbrot(ca,cb);
            v = Math.floor(v*(255/maxIteracij));
            ctx.fillStyle = "rgb("+v+","+v+","+v+")";
            ctx.fillRect(x,y,1,1);
        }
    }
}

Da bodo piksli, ki so v mandelbrotovi množici, pobarvani z belo, vrednost funkcije pomnožimo z izrazom 255/maxIteracij. Pri tem je nujno, da vrednost izraza zaokrožimo, ker format barve rgb() ne sprejema vrednosti z decimalkami.

Če zdaj zaženemo program in pokličemo funkcijo izris(), dobimo črn zaslon z majhno belo piko (na sliki je povečana, da se jo vidi, v resnici je velika le 1px) v zgornjem levem kotu. To se zgodi zato, ker je izhodišče našega koordinatnega sistema v zgornjem levem kotu in se razteza po \(x\) osi do 600 in po \(y\) osi do 400, mandelbrotova množica pa se nahaja med \(c_a\) od -2 do 1 in \(c_b\) od -1 do 1.

Da to popravimo moramo \(c_a\) in \(c_b\) prilagoditi, da bosta pribljižno v intervalu od -2 do 2. Za začetek lahko delimo x in y z vrednostjo zoom (recimo, da je zoom=100), s tem dosežemo, da se novi x rasteza med 0 in 6 in y med 0 in 4.

\(x∈(0,600)\)                               \(y∈(0,400)\)

\(c_a=x/k_{zoom}\)      \(c_a∈(0,6)\)      \(c_b=y/k_{zoom}\)      \(c_b∈(0,4)\)

Če bi zdaj zagnali program, bi videli en del mandelbrotove množice v zgornjem levem kotu, vendar bi bil velik del odrezan. Da to popravimo moramo koordinatno izhodišče prestaviti na sredino zaslona. To naredimo tako, da od x odštejemo polovico širine zaslona, od y pa polovico višine zaslona.

\(c_a=\frac{x-w/2}{k_{zoom}}\)      \(c_a∈(-3,3)\)      \(c_b=\frac{y-h/2}{k_{zoom}}\)      \(c_b∈(-2,2)\)

Zdaj \(c_a\) in \(c_b\) tesneje zajemata mandelbrotovo množico. Ampak kaj pa če bi želeli koordinatno izhodišče (oziroma sliko fraktala) premakniti? To naredimo tako, da od x odštejemo zamikX in od y odštejemo zamikY. Če je \(x_{zamik}=100\) in \(y_{zamik}=50\), dobimo intervala:

\(c_a=\frac{x-w/2-x_{zamik}}{k_{zoom}}\)      \(c_a∈(-4,2)\)      \(c_b=\frac{y-h/2-y_{zamik}}{k_{zoom}}\)      \(c_b∈(-2.5,1.5)\)

To bo uporabno, ko bomo želeli povečati fraktal in se premikati po njem. Zdaj lahko enačbo vstavimo v program:

function izris(){
    for(var x=0; x<w; x++){
        for(var y=0; y<h; y++){
            var ca = (x-w/2-zamikX)/zoom;
            var cb = (y-h/2-zamikY)/zoom;
            var v = mandelbrot(ca,cb);
            v = Math.floor(v*(255/maxIteracij));
            ctx.fillStyle = "rgb("+v+","+v+","+v+")";
            ctx.fillRect(x,y,1,1);
        }
    }
}

Zdaj lahko fraktal povečamo in premaknemo na sredino.

Dodajanje kontrol za raziskovanje fraktala

Da se bomo lahko premikali po fraktalu, in ga povečali dodajmo nekaj kode, ki bo ob pritisku tipke spremenila pozicijo ali povečavo fraktala.

window.addEventListener('keydown', function(event){
    if(event.code == "ArrowUp"){
        zamikY += 100;
    }
    else if(event.key == "ArrowDown"){
        zamikY -= 100;
    }
    else if(event.key == "ArrowLeft"){
        zamikX += 100;
    }
    else if(event.key == "ArrowRight"){
        zamikX -= 100;
    }
    else if(event.key == "+"){
        zoom *= 2;
        zamikX *= 2;
        zamikY *= 2;
    }
    else if(event.keyCode == "-"){
        zoom /= 2;
        zamikX /= 2;
        zamikY /= 2;
    }
    izris();
});

Da se nam fraktal ne odmakne, ko ga povečujemo/pomanjšujemo moramo pri povečevanju/pomanjševanju povečati/zmanjšati tudi zamik.

S tem smo program končali, jaz ga bom še dopolnil (barvanje fraktala, prilagajanje velikosti okna, nastavljanje natančnosti …) ampak tega ne bom naprej razlagal, če vas zanima končna verzija programa si lahko ogledate izvorno kodo ali pa preizkusite demo.

preizkusi program

Za za premikanje uporabljajte puščice na tipkovnici, za zoom pa tipki + in -.

Rezultati

Tukaj je nekaj slik, ki sem jih naredil z opisanim programom. Način barvanja sem spremenil iz sivinske v barvno paleto (črna, rdeča, zelena, modra, črna).

Omejitve programa

Največja omejitev tega programa je hitrost. Da program za vsak piksel ugotovi ali je v mandelbrotovi množici, traja zelo dolgo časa. Ta problem lahko rešimo tako, da zmanjšamo število pikslov ali pa zmanjšamo število iteracij, v obeh primerih izgubimo na natančnosti. Druga omejitev je natančnost računanja. JavaScript uporablja 64-bit števila s plavajočo vejico, kar pomeni, da bo pri računanju zelo majhnih števil naredil napako. (vir: https://www.w3schools.com/js/js_numbers.asp) Zato bomo videli nepravilnosti, ko bomo zelo povečali fraktal.

Na sliki so prikazane nepravilnosti, ki se pojavijo, ko fraktal zelo povečamo.

Pohitritev programa z WebGL

WebGl je API, ki omogoča da z programskim jezikom JavaScript, uporabljamo grafično kartico. Grafični procesor ima veliko več jeder, kot CPE zato je računanje mandelbrotove množice veliko hitrejše, tudi pri višji ločljivosti. Žal ima WebGL (vsaj ternutan verzija) hudo pomankljivost. Podpira samo 32 bitna števila, kar pomeni da še hitreje pride do nepravilnosti.

Tukaj lahko preizkusite program ki uporablja WebGL:

preizkusi program

Za premikanje po fraktalu lahko uporabljate tipke na tipkovnici, miško ali pa zaslon na dotik.

Datoteke