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).
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.
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.
Obrázek 3: Logistická mapaPří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.
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.
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;
}
}
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í.
Obrázek 7: První ukázka fraktálu vytvořeného ve druhém demonstračním příkladu Obrázek 8: Druhá ukázka fraktálu vytvořeného ve druhém demonstračním příkladu Obrázek 9: Třetí ukázka fraktálu vytvořeného ve druhém demonstračním příkladu 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
- Stránky programu FractInt:
http://www.fractint.org - Aleksandr Lyapunov (Wikipedia):
http://en.wikipedia.org/wiki/Aleksandr_Lyapunov - Lyapunov exponent (Wikipedia):
http://en.wikipedia.org/wiki/Lyapunov_exponent - Lyapunov fractal (Wikipedia):
http://en.wikipedia.org/wiki/Lyapunov_fractal - The Lyapunov exponent:
http://monet.unibas.ch/~elmer/pendulum/lyapexp.htm - Lyapunov Fractals:
http://spanky.triumf.ca/www/fractint/lyapunov_type.html - Lyapunov Space:
http://hypertextbook.com/chaos/44.shtml - Markus-Lyapunov Fractals:
http://perso.orange.fr/charles.vassallo/en/lyap_art/lyapdoc.html - Logistic function:
http://en.wikipedia.org/wiki/Logistic_function - Logistic map:
http://en.wikipedia.org/wiki/Logistic_map - The Logistic Map:
http://www.geocities.com/CapeCanaveral/Hangar/7959/logisticmap.html - List of chaotic maps:
http://en.wikipedia.org/wiki/List_of_chaotic_maps
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ů.
Obrázek 12: Detail jiného Markus-Ljapunovova fraktálu