Hlavní navigace

Markus-Ljapunovovy fraktály

Pavel Tišnovský

Dnes si ukážeme tvorbu Markus-Ljapunovových fraktálů, pomocí kterých je možné vytvářet zajímavé obrázky použitelné například jako textury či výšková pole. Algoritmus pro výpočet těchto fraktálů je možné různě modifikovat; my si na příkladech ukážeme dvě rozdílné verze, na které později navážeme dalším rozšířením.

Obsah

1. Markus-Ljapunovovy fraktály
2. Bifurkační diagram a logistická mapa
3. Numerický výpočet Ljapunovova exponentu a Markus-Ljapunovův fraktál
4. První demonstrační příklad – základní Markus-Ljapunovův fraktál vykreslený ve stupních šedi
5. Markusovo a Hessovo rozšíření v podobě výběru posloupnosti funkcí
6. Druhý demonstrační příklad – Markusovo a Hessovo rozšíření spolu s aplikací barvové palety
7. Odkazy na další informační zdroje
8. Obsah dalšího pokračování tohoto seriálu

1. Markus-Ljapunovovy fraktály

Markus-Ljapunovovy fraktály tvoří poslední skupinu fraktálů vykreslovaných v rovině, kterou se budeme v tomto seriálu zabývat. Jedná se o fraktály vytvářené v pravidelné rastrové mřížce většinou zobrazované ve formě bitmapy či pixmapy (podobně jako známá Mandelbrotova množina a množiny Juliovy), ve které je každému pixelu přiřazena barva v závislosti na numericky vypočtené hodnotě takzvaného Ljapunovova exponentu. Pomocí tohoto exponentu je v matematice i fyzice určována (velmi zjednodušeně řečeno) míra lineární stability či naopak nestability atraktoru nějakého dynamického systému pro dané vstupní podmínky; dokonce se jedná o jeden z nejdůležitějších parametrů mnoha nelineárních dynamických systémů (včetně modelů počasí, modelů populačního růstu, který si blíže popíšeme, modelů pohybů cen na burze atd).

fractals69_1
Obrázek 1: Ukázka Markus-Ljapunovova fraktálu vytvořeného v programu FractInt

Ljapunovův exponent je poměrně složité analyticky vypočítat (je počítán pomocí limity pro t jdoucí k nekonečnu), v dále uváděných praktických příkladech si však vystačíme s prostým numerickým přiblížením přesného výpočtu. Markus-Ljapunovovy fraktály jsou odvozeny z jednorozměrných bifurkačních diagramů a logistické mapy – oba dva termíny jsme si (spolu s několika praktickými příklady) popsali v prvních částech tohoto seriálu; ve druhé kapitole si některé základní informace zopakujeme. Zbývá dodat důležitou poznámku (zejména pro vyhledávání dokumentů pojednávajících o Ljapunovových exponentech a fraktálech na Internetu): Ljapunov je slovanské jméno původně psané azbukou ляпунов, které se většinou v anglické transkripci přepisuje jako Lyapunov.

fractals69_2
Obrázek 2: Další ukázka téhož typu fraktálu vytvořeného v programu FractInt

2. Bifurkační diagram a logistická mapa

Při studiu populačního růstu v uzavřeném prostoru (například se může jednat o ryby žijící v jednom rybníce) se pro jeho značně zjednodušený model došlo k závěru, že velikost populace v určitém časovém období, například pro každý rok, závisí na velikosti populace v jednom předchozím časovém období. Systém tedy nemá žádnou dlouhodobou „paměť“, ve které by si pamatoval všechny předešlé stavy. Populační růst klesá, jakmile celková část populace dosáhne či přesáhne určitou hodnotu Xm (příliš mnoho organismů v uzavřeném prostoru), a v menší populaci naopak dochází k jejímu růstu. To je způsobeno tím, že populace by sice principielně rostla geometrickou řadou, ale zdroje nutné k přežití (potrava, životní prostor atd.) klesají kvadraticky, proto se původní neomezený růst populace zastavuje a dokonce může dojít k jejímu poklesu.

Dynamický systém takto navrženého modelu populace (dynamický je proto, že se mění v čase) je nelineární a současně chaotický, vžil se pro něj také název Verhulstův proces. Ve skutečnosti se jedná o velmi jednoduchý nelineární dynamický systém popsaný jednou rovnicí, v praxi se můžeme setkat se systémy mnohem složitějšími (chemické reakce, biologické procesy, fyzikální jevy atd). Výpočet populace v n+1 generaci vychází z velmi jednoduchého matematického modelu omezeného růstu populace po jednotlivých generacích:

xn+1=Gr × xn × (1-xn)

Ve výše uvedeném vzorci se symbolem xn značí počet jedinců v n-té generaci, symbolem xn+1 počet jedinců v následující generaci a Gr je velikost růstu (resp. míra růstu). Hodnoty xn se zadávají v rozmezí 0..1, kde 0 značí vymření populace a 1 plně saturovaný stav s ohledem na dostupné zdroje.

V pátém pokračování tohoto seriálu jsme si také řekli, že pro některé hodnoty růstu Gr se populace po několika krocích ustálí na nějaké konstantní hodnotě; jiné hodnoty zapříčiní oscilaci mezi dvěma, čtyřmi, osmi atd. stavy a konečně pro většinu hodnot Gr větších než 3,57 se systém stává chaotickým, tj. není možné předpovědět jeho budoucí stav (pomocí podmínky pro stacionární body, tj. body, pro které platí xn+1=xn, lze jednoduše zjistit, že systém je vždy stabilní pro Gr ležící v rozsahu –1..3). I pro některé hodnoty růstu větší než „magické číslo“ 3,57 však můžeme nacházet oblasti, ve kterých je systém stabilní tj. po ustálení periodický, vždy s periodami o hodnotě mocniny dvou. Přímým zobrazením hodnot xn získáme takzvanou logistickou mapu.

fractals69_3
Obrázek 3: Logistická mapa

Přímé zobrazení logistické funkce nám však nedává zcela jasnou představu o tom, při jakých počátečních podmínkách je systém stabilní, kdy dochází ke zdvojování period (bifurkaci) a pro jaké hodnoty nastává chaotické chování systému. Pro jasnější představu o chování jednorozměrného dynamického systému se proto používá jiný účinný vizuální nástroj – bifurkační diagram. V bifurkačním diagramu jsou na horizontální osu naneseny hodnoty některého vstupního parametru (v našem případě růstu Gr) a na vertikální osu vnitřní stavy dynamického systému, což je u logistické mapy hodnota xn, jichž systém nabývá pro daný počet kroků (iterací).

Iterovaná hodnota xn je vykreslena triviálním způsobem pomocí jednoho bodu, nedochází tedy k tvorbě spojitého grafu (vzhledem k jeho charakteru to ani není možné). V některých případech se bifurkační diagram dále upravuje, například tak, že se vykreslují hodnoty až po předem zadaném počtu iterací, kdy je systém již částečně ustálen (jedná se o takzvané „počáteční“ nebo „startovní“ iterace) a bifurkační diagram je v tomto případě označen jako filtrovaný bifurkační diagram.

fractals69_4
Obrázek 4: Bifurkační diagram

3. Numerický výpočet Ljapunovova exponentu a Markus-Ljapunovův fraktál

Vzhledem k tomu, že se Ljapunovův exponent vyjadřuje analyticky, je poměrně obtížné ho v této přímé podobě programově vypočítat. Můžeme si však pomoci (obecně méně přesným) numerickým výpočtem, který je pro programovou implementaci jednoduchý a přímočarý. Ukažme si to na příkladu výše popisovaného modelu růstu definovaného vztahem:

xn+1=Gr × xn × (1-xn)

Pro tento vzorec by se Ljapunovův exponent vyjádřil jako suma hodnot pro všechna xn:

když |Gr-2Grxn|>0 pak:
suma=suma+log2|Gr-2Grxn|

Pravděpodobně jste si všimli, že výraz v absolutní hodnotě odpovídá derivaci pravé strany původního vzorce podle xn. Pro dostatečně velký počet průchodů (v každém průchodu se zvýší index n) dostaneme příslušnou hodnotu Ljapunovova exponentu. Pokud je tento exponent záporný, jedná se o stabilní oblast, pokud je kladný, jde o oblast s chaotickým chováním.

Markus-Ljapunovův fraktál přímo vychází z výpočtu Ljapunovova exponentu, pouze ho nápaditým způsobem rozšiřuje do dvou dimenzí. Vzhledem k tomu, že se má vytvořit plošný obrázek (bitmapa, pixmapa), použijí se při výpočtech dva dynamické systémy modelu růstu, přičemž se v každém pixelu zkoumá Ljapunovův exponent pro jinou hodnotu míry růstu Gr. Ve směru horizontální osy se zvětšuje míra růstu prvního dynamického systému, ve směru osy vertikální pak míra růstu druhého systému. V demonstračních příkladech jsou tyto hodnoty uloženy v proměnných nazvaných startx a starty, samotná hodnota xn je pak uložena v proměnné z, aby se zdůraznilo, že se vlastně jedná o třetí souřadnici odpovídající barvě pixelu.

Při iterativním výpočtu Ljapunovova exponentu se v iterační smyčce pravidelně střídají vzorce obou dynamických systémů, což v důsledku vede k tvorbě zdánlivě „trojrozměrných“ obrázků. V případě, že je výsledná hodnota Ljapunovova exponentu kladná (chaos), je daný pixel vykreslený konstantní barvou (například černou), v opačném případě je obarven na základě vypočtené hodnoty exponentu. Do celého výpočtu je ještě vložena podmínka, která zabrání růstu xn nade všechny meze. Vše si prakticky ukážeme v následující kapitole.

fractals69_5
Obrázek 5: Markus-Ljapunovův fraktál vykreslený v programu FractInt

4. První demonstrační příklad – základní Markus-Ljapunovův fraktál vykreslený ve stupních šedi

Dnešní první demonstrační příklad slouží k výpočtu Markus-Ljapunovova fraktálu spolu s jeho vykreslením do pixmapy, ve které jsou jednotlivé pixely obarveny pouze různými úrovněmi šedé barvy (při výpočtech tedy není aplikováno žádné mapování pomocí barvové palety). Každému pixelu ve výsledné pixmapě je přiřazena dvojice počátečních hodnot označených startx a starty, které odpovídají velikostem růstu (v původním vzorci označeném jako Gr). Tyto hodnoty leží v rozsahu <-8, 8> pro startx a <-6, 6> pro starty.

Při výpočtu Ljapunovova exponentu se provádí iterační výpočet s proměnnou z. V tomto výpočtu se pravidelně střídají dvě logistické funkce, jedna s růstem startx a druhá s růstem starty. Zdrojový text prvního demonstračního příkladu můžete získat zde, popř. i jako HTML soubor se zvýrazněnou syntaxí. Funkce provádějící výpočet barev pixelů v tomto fraktálu má tvar:

//-----------------------------------------------------------------------------
// Vypocet Markus-Ljapunovovych fraktalu
//-----------------------------------------------------------------------------
void recalcLyapunov(pixmap *pix,                // pracovni pixmapa
                    int    maxiter,             // maximalni pocet iteraci
                    double xmin,                // meze fraktalu v rovine
                    double ymin,
                    double xmax,
                    double ymax)
{
    int x, y, i, j;                             // pocitadla smycek
    double startx, starty;
    double z, r=0.0;
    double total;

    starty=ymin;
    // pro vsechny radky pixmapy
    for (y=0; y<pix->height; y++) {
        startx=xmin;
        printf("%d\t", y);
        // pro vsechny pixely na radku
        for (x=0; x<pix->width; x++) {
            startx+=(xmax-xmin)/pix->width;     // pocatecni hodnota prvni funkce
            z=0.5;
            total=0.0;
            j=0;                                // pocitadlo pro pristup k retezci
            // iteracni smycka
            for (i=0; i<maxiter; i++) {
                switch (i%2) {                  // vyber iteracni funkce
                    case 0:
                        z=(starty*z)*(1.0-z);
                        r=starty;
                        break;
                    case 1:
                        z=(startx*z)*(1.0-z);
                        r=startx;
                        break;
                }
                j++;
                if (i>maxiter/2) {              // ve druhe polovine cyklu se pocita Ljapunuv exponent
                    if (fabs(r-2.0*r*z)>0.0)
                        total+=log(fabs(r-2.0*r*z))/log(2.0);
                }
                if (z>500.0) break;
            }
            total=150.0*total/maxiter;          // vypocet exponentu s vynasobenim konstantou 150
            if (total<0) {                      // obarvit pouze pixely se zapornym exponentem
                if (total<-255) total=-255;
                total=-total/255;               // dalsi tri programove radky - gamma korekce
                total=pow(total, 0.5);
                total*=255.0;
                putpixel(pix, x, y, total, total, total);
            }
            else                                // "chaoticke" pixely jsou cerne
                putpixel(pix, x, y, 0, 0, 0);
        }
        // prechod na dalsi radek v pixmape
        starty+=(ymax-ymin)/pix->height;
    }
} 

fractals69_6
Obrázek 6: Screenshot prvního demonstračního příkladu

5. Markusovo a Hessovo rozšíření v podobě výběru posloupnosti funkcí

Markus a Hess rozšířili výše popsaný výpočet Markus-Ljapunovova fraktálu o další možnost, resp. o další stupeň volnosti. Místo pevně dané posloupnosti funkcí růstu při iterativním výpočtu Ljapunovova exponentu, které můžeme označit například symboly a a b, je možné použít posloupnost libovolnou. Specifikace vybrané posloupnosti může být zadána buď pomocí řetězce (například řetězec „abba“ značí postupný výběr funkcí a, b, b, a, a, b, b, a …) nebo i decimálním číslem, které vznikne po převodu binárních symbolů (0=a, 1=a) do sekvence bitů chápané jako binární zápis celého nezáporného čísla.

Tímto způsobem je zápis posloupnosti řešen například v programu FractInt, konkrétně se jedná o parametr Order u typu fraktálu nazvaného Lyapunov. Výběr funkce z posloupnosti je jednoduchý – při iterativním výpočtu se použije index ukazující na znak v řetězci, kterým se daná funkce volí. Původní posloupnost (použitá mimochodem i v dnešním prvním demonstračním příkladu) by se zapsala řetězcem o dvou znacích „ab“, protože se při dosažení konce řetězce automaticky přejde zpět na jeho počátek. Céčkovská funkce, která umožňuje vypočítat a nakreslit rozšířený Markus-Ljapunovův fraktál, může (při vynechání různých pro ukázku nepodstatných optimalizačních postupů) vypadat následovně:

//-----------------------------------------------------------------------------
// Vypocet Markus-Lyapunovych fraktalu
//-----------------------------------------------------------------------------
void recalcLyapunov(pixmap *pix,                // pracovni pixmapa
                    int    palette,             // cislo barvove palety
                    int    maxiter,             // maximalni pocet iteraci
                    double xmin,                // meze fraktalu v rovine
                    double ymin,
                    double xmax,
                    double ymax,
                    const char *order)          // poradi aplikace funkci
{
    int x, y, i, j;                             // pocitadla smycek
    double startx, starty;
    double z, r=0.0;
    double total;

    starty=ymin;
    // pro vsechny radky pixmapy
    for (y=0; y<pix->height; y++) {
        startx=xmin;
        printf("%d\t", y);
        // pro vsechny pixely na radku
        for (x=0; x<pix->width; x++) {
            startx+=(xmax-xmin)/pix->width;     // pocatecni hodnota prvni funkce
            z=0.5;
            total=0.0;
            j=0;                                // pocitadlo pro pristup k retezci
            // iteracni smycka
            for (i=0; i<maxiter; i++) {
                switch (order[j]) {             // vyber iteracni funkce
                    case 'a':
                        z=(starty*z)*(1.0-z);
                        r=starty;
                        break;
                    case 'b':
                        z=(startx*z)*(1.0-z);
                        r=startx;
                        break;
                    default:                    // spatny znak v retezci
                        z=z;
                        r=0;
                }
                j++;
                if (j>strlen(order)) j=0;       // postupne projit celym retezcem order

                if (i>maxiter/2) {              // ve druhe polovine cyklu se pocita Lyapunuv exponent
                    if (fabs(r-2.0*r*z)>0.0)
                        total+=log(fabs(r-2.0*r*z))/log(2.0);
                }
                if (z>500.0) break;
            }
            total=255.0*total/maxiter;          // uprava ukazatele do barvove palety
            if (total<0) {                      // obarvit pouze pixely se zapornym exponentem
                unsigned char r, g, b;
                if (total<-255) total=-255;
                mapPalette(palette, -total, &r, &g, &b); // vyber barvy z barvove palety
                putpixel(pix, x, y, r, g, b);   // vykresleni pixelu
            }
            else                                // pro "chaoticke" posloupnosti jsou pixely cerne
                putpixel(pix, x, y, 0, 0, 0);
        }
        // prechod na dalsi radek v pixmape
        starty+=(ymax-ymin)/pix->height;
    }

} 

6. Druhý demonstrační příklad – Markusovo a Hessovo rozšíření spolu s aplikací barvové palety

V předchozí kapitole jsme si uvedli céčkovskou funkci, kterou je možné použít pro vykreslování Hessem rozšířených Markus-Ljapunovových fraktálů. Jako jeden z parametrů této funkce se předává řetězec, ve kterém je pomocí znaků a a b zapsána posloupnost funkcí růstu tak, jak jsou tyto funkce aplikovány při iterativním výpočtu. Ostatní znaky v tomto řetězci se ignorují. Posloupnost by měla obsahovat alespoň jeden znak a a jeden znak b, v opačném případě dostaneme nezajímavý „jednorozměrný“ obrázek obsahující pouze gradientní barvové přechody orientované v horizontálním či vertikálním směru.

Příklad je také rozšířen o možnost výběru barvové palety, pomocí které jsou počítané pixely obarveny. Na obrázcích 7–10 jsou ukázány příklady fraktálů vykreslených pomocí druhého demonstračního příkladu. Zdrojový text druhého demonstračního příkladu je možné získat zde, popř. i jako HTML soubor se zvýrazněnou syntaxí.

fractals69_7
Obrázek 7: První ukázka fraktálu vytvořeného ve druhém demonstračním příkladu

fractals69_8
Obrázek 8: Druhá ukázka fraktálu vytvořeného ve druhém demonstračním příkladu

fractals69_9
Obrázek 9: Třetí ukázka fraktálu vytvořeného ve druhém demonstračním příkladu

fractals69_a
Obrázek 10: Čtvrtá ukázka fraktálu vytvořeného ve druhém demonstračním příkladu

7. Odkazy na další informační zdroje

  1. Stránky programu FractInt:
    http://www.frac­tint.org
  2. Aleksandr Lyapunov (Wikipedia):
    http://en.wiki­pedia.org/wiki/A­leksandr_Lyapu­nov
  3. Lyapunov exponent (Wikipedia):
    http://en.wiki­pedia.org/wiki/Ly­apunov_exponent
  4. Lyapunov fractal (Wikipedia):
    http://en.wiki­pedia.org/wiki/Ly­apunov_fractal
  5. The Lyapunov exponent:
    http://monet.u­nibas.ch/~elmer/pen­dulum/lyapexp­.htm
  6. Lyapunov Fractals:
    http://spanky­.triumf.ca/www/frac­tint/lyapunov_ty­pe.html
  7. Lyapunov Space:
    http://hypertex­tbook.com/cha­os/44.shtml
  8. Markus-Lyapunov Fractals:
    http://perso.o­range.fr/char­les.vassallo/en/ly­ap_art/lyapdoc­.html
  9. Logistic function:
    http://en.wiki­pedia.org/wiki/Lo­gistic_functi­on
  10. Logistic map:
    http://en.wiki­pedia.org/wiki/Lo­gistic_map
  11. The Logistic Map:
    http://www.ge­ocities.com/Ca­peCanaveral/Han­gar/7959/logis­ticmap.html
  12. List of chaotic maps:
    http://en.wiki­pedia.org/wiki/Lis­t_of_chaotic_maps

fractals69_b
Obrázek 11: Detail Markus-Ljapunovova fraktálu

8. Obsah dalšího pokračování tohoto seriálu

V následucím pokračování tohoto seriálu dokončíme část věnovanou Markus-Ljapunovovým fraktálům. Ukážeme si rozšíření těchto fraktálů o další funkce, které mohou nahradit původní funkci populačního růstu a pomocí nichž je možné zajistit vytváření mnoha dalších zajímavých obrázků.

fractals69_c
Obrázek 12: Detail jiného Markus-Ljapunovova fraktálu
Našli jste v článku chybu?

28. 2. 2007 10:35

Diky za opravu, azbuku sice porad znam (presneji, naposled jsem ji aktivne pouzival na zakladce), ale uz jsem se pri hledani Unikodu zapomel podivat na velikost - holt tam nejsou tak viditelne rozdily jako v latince.

28. 2. 2007 9:13

Honza (neregistrovaný)
V azbuce se taky používá velké Л, takže Ляпунов. Celým jménem Александр Михайлович Ляпунов (Alexandr Michajlovič L., 1857-1918).
Vitalia.cz: Jmenuje se Janina a žije bez cukru

Jmenuje se Janina a žije bez cukru

Měšec.cz: Kdy vám stát dá na stěhování 50 000 Kč?

Kdy vám stát dá na stěhování 50 000 Kč?

DigiZone.cz: Recenze Westworld: zavraždit a...

Recenze Westworld: zavraždit a...

Lupa.cz: Teletext je „internetem hipsterů“

Teletext je „internetem hipsterů“

Podnikatel.cz: Změny v cestovních náhradách 2017

Změny v cestovních náhradách 2017

Měšec.cz: Jak vymáhat výživné zadarmo?

Jak vymáhat výživné zadarmo?

Lupa.cz: Propustili je z Avastu, už po nich sahá ESET

Propustili je z Avastu, už po nich sahá ESET

Podnikatel.cz: Udávání kvůli EET začalo

Udávání kvůli EET začalo

Root.cz: Certifikáty zadarmo jsou horší než za peníze?

Certifikáty zadarmo jsou horší než za peníze?

Měšec.cz: U levneELEKTRO.cz už reklamaci nevyřídíte

U levneELEKTRO.cz už reklamaci nevyřídíte

Podnikatel.cz: Chaos u EET pokračuje. Jsou tu další návrhy

Chaos u EET pokračuje. Jsou tu další návrhy

Lupa.cz: Insolvenční řízení kvůli cookies? Vítejte v ČR

Insolvenční řízení kvůli cookies? Vítejte v ČR

Podnikatel.cz: Chtějte údaje k dani z nemovitostí do mailu

Chtějte údaje k dani z nemovitostí do mailu

Vitalia.cz: Pamlsková vyhláška bude platit jen na základkách

Pamlsková vyhláška bude platit jen na základkách

Vitalia.cz: Proč vás každý zubař posílá na dentální hygienu

Proč vás každý zubař posílá na dentální hygienu

Vitalia.cz: Potvrzeno: Pobyt v lese je skvělý na imunitu

Potvrzeno: Pobyt v lese je skvělý na imunitu

Podnikatel.cz: Podnikatelům dorazí varování od BSA

Podnikatelům dorazí varování od BSA

Měšec.cz: Zdravotní a sociální pojištění 2017: Připlatíte

Zdravotní a sociální pojištění 2017: Připlatíte

120na80.cz: Co všechno ovlivňuje ženskou plodnost?

Co všechno ovlivňuje ženskou plodnost?

Lupa.cz: Co se dá měřit přes Internet věcí

Co se dá měřit přes Internet věcí