Hlavní navigace

Fraktály v počítačové grafice XXII

Pavel Tišnovský 21. 3. 2006

V dnešním pokračování seriálu o fraktálech používaných (nejenom) v počítačové grafice si popíšeme další velmi zajímavý fraktál, který je možné vytvořit v komplexní rovině. Jedná se o takzvanou Newtonovu fraktální množinu, která vzniká barevným zvýrazněním obecných komplexních kořenů jednoduchých polynomů.

Obsah

1. Newtonova iterační metoda pro obor reálných čísel
2. Nalezení kořenů komplexního polynomu
3. Citlivost na počáteční podmínky pro polynomy třetího a vyššího stupně
4. Odvození funkce Newtonovy iterační metody pro polynom z3-1=0
5. Newtonův fraktál polynomu z3-1=0
6. Modifikovaný Newtonův fraktál polynomu z3-1=0
7. Newtonovy fraktály polynomů vyšších stupňů
8. Obsah dalších pokračování tohoto seriálu

1. Newtonova iterační metoda pro obor reálných čísel

Jedním z velmi častých úkolů, se kterými se setkáváme v mnoha vědních disciplínách, je nalezení kořenů polynomů, tj. takových hodnot x (reálných či komplexních), při jejichž zpětném dosazení do rovnice polynomu anxn+an-1xn-1…a2x2+a1x+a0 vyjde nulový výsledek – to znamená, že graf polynomické funkce v tomto bodě protíná x-ovou osu. Pro hledání kořenů polynomů je možné použít značně velké množství numerických metod, které se liší svou přesností, stabilitou a mimo jiné také rychlostí vyhledání kořenů, protože absolutně přesné analytické metody jsou odvozeny pouze pro malé hodnoty n.

Jednou z nejpoužívanějších a současně i nejznámějších numerických metod určených pro vyhledávání kořenů polynomů je Newtonova iterační metoda, která se stala populární zejména díky tomu, že se numerické řešení poměrně rychle přibližuje ke kořenu, v mnoha případech je pro dosažení požadované přesnosti nutné provést pouze několik iterací, přičemž v každé iteraci se výsledek zpřesňuje o jedno až dvě desetinná místa. Newtonova iterační metoda, kterou je možné jednoduše vysvětlit i geometricky (ve své podstatě vypočítává průsečík tečny v nějakém bodě grafu funkce s x-ovou osou), spočívá v opakovaném výpočtu rovnice:

xi+1=xi-f(xi)/f'(xi)

přičemž je levá strana rovnice (ve výše uvedeném vzorci označena jako xi+1) po výpočtu opět dosazena na vstup. Newtonova iterační metoda pro daný polynom vždy nalezne pouze jeden kořen polynomu. To, jaký kořen je nalezen, závisí především na počáteční hodnotě x0, protože je nalezen nejbližší kořen k této hodnotě. Jedná se také o metodu, kterou je možné použít nejenom u polynomů, ale i u všech dalších funkcí, které mají v každém bodě derivaci, přičemž by derivace neměla být nulová – to může nastat v lokálních maximech a minimech, ve kterých se řešení pomocí této iterační metody „zasekne“ (nutno však dodat, že i další metody založené na postupné iteraci trpí podobným problémem).

Jak se můžeme jednoduchými geometrickými i numerickými pokusy přesvědčit, je nalezení kořenů polynomů v oboru reálných čísel poměrně přímočaré a nevyskytují se zde žádné neočekávané nepravidelnosti (geometrické pokusy se dají provádět na zobrazeném grafu funkce s pomocí pravítka, kterým se odhadne směrnice a nalezne se průsečík směrnice s x-ovou osou grafu). Reálná resp. x-ová osa je rozdělena do několika na sebe navazujících oblastí, kde každá oblast odpovídá jednomu kořenu. Pokud se počáteční hodnota x nachází v některé z těchto oblastí, řešení rychle dospěje k odpovídajícímu kořenu. Pojďme se však nyní podívat, jak celá situace vypadá v komplexní rovině.

fractals21_1
Obrázek 1: Princip Newtonovy iterační metody

2. Nalezení kořenů komplexního polynomu

Mnozí ze čtenářů se pravděpodobně podivují, proč je v seriálu o fraktálech popisována Newtonova iterační metoda, ta přece nemá (alespoň na první pohled) s fraktály nic společného. Ukažme si však nyní použití Newtonovy iterační metody v komplexní rovině, ze kterého vztah k fraktálům za malou chvíli vyplyne. Platí, že každý polynom stupně n má také n kořenů, přičemž tyto kořeny mohou být i násobné, tj. více kořenů může mít stejnou hodnotu.

Kořeny polynomu mohou být v některých případech reálné, tj. budou ležet na reálné ose. Reálné kořeny má například polynom x2-1=0. U velké části polynomů však leží některé nebo dokonce všechny kořeny v komplexní rovině. Typickým příkladem může být polynom x2+1=0, který sice nemá v oboru reálných čísel žádný kořen (rovnici nelze v tomto oboru řešit), ale v komplexní rovině je samozřejmě možné kořeny nalézt, zejména když si uvědomíme platnost známého vztahu i2=-1. Newtonovu iterační metodu je možné použít i v oboru komplexních čísel, příslušný iterační vztah se přitom změní pouze nepatrně:

zi+1=zi-f(zi)/f'(zi)

Mohlo by se zdát, že všechny výše uvedené vlastnosti Newtonovy iterační metody zůstanou zachovány, například i vlastnost nalezení nejbližšího kořenu k zadané hodnotě z0. Ve skutečnosti to však není zdaleka pravda, o čemž se můžeme přesvědčit v následující kapitole.

3. Citlivost na počáteční podmínky pro polynomy třetího a vyššího stupně

Pro polynom z3-1=0 je možné analyticky spočítat všechny tři kořeny: jsou to komplexní hodnoty 1+i0, cos 2π/3 + i sin 2π/3 a cos 4π/3 + i sin 4π/3. Pokud budeme do komplexní roviny vykreslovat body z0 obarvené jednou ze tří barev podle toho, ke kterému kořenu řešení dospěje, získáme vlastně mapu atraktorů daného polynomu – pokud se bod z0 nachází v ploše obarvené jednou barvou, bude odpovídající řešení „přitahováno“ ke kořenu této barvy. Mohli bychom se domnívat, že celá komplexní rovina bude rozdělena na tři stejně velké barevné oblasti s pravidelnými hranicemi; jednalo by se o nekonečně velké výseče střetávající se v počátku.

Skutečnost je však mnohem složitější (a zajímavější), protože při obarvování bodů dospějeme k fraktálnímu tvaru, který je zobrazen na druhém obrázku. Všimněte si, že komplexní rovina je skutečně rozdělena na tři velké různobarevné plochy, kde každá plocha odpovídá jednomu řešení/kořenu. Mezi těmito plochami se však nenachází pravidelná hranice, ale naopak fraktálovitě se členící a soběpodobná oblast, kde dochází k nekonečnému dělení na podoblasti. Detailnější pohled na rekurzivní dělení je zobrazen na třetím obrázku.

fractals22_2
Obrázek 2: Barevně vyznačené oblasti působnosti tří atraktorů

Fraktální charakteristika hranice mezi třemi velkými barevnými oblastmi mimo jiné naznačuje značnou citlivost Newtonovy metody použité pro komplexní polynomy – dva nekonečně blízké počáteční body z0A a z0B se mohou po několika iteracích přiblížit ke dvěma rozdílným atraktorům (kořenům).

Přitom je zajímavé, že podobný obrázek získáme při vizualizaci výsledků měření jednoduchého pokusu se třemi magnety a zavěšeným závažím: magnety jsou souhlasným pólem připevněny na podložku tak, že tvoří vrcholy rovnostranného trojúhelníku. Přesně nad střed tohoto trojúhelníku je přivázáno železné či ocelové závaží (provázek by měl být co nejdelší, aby se eliminovala gravitace). Rukou můžeme závaží vychýlit z rovnovážné polohy. Po opětovném puštění závaží dojde k jeho ustálení nad jedním z magnetů, samozřejmě za předpokladu, že jsou magnety dostatečně silné, aby jejich přitažlivá síla překonala gravitační sílu, která se závaží snaží ponechat ve svislé poloze.

Mohli bychom očekávat, že pokud budeme obarvovat místa, ze kterých bylo závaží puštěno, podle toho, ke kterému magnetu se nakonec vychýlí, vzniknou tři pravidelné barevné oblasti. Ve skutečnosti je však systém závaží–magnety velmi citlivý na počáteční podmínky a na hranicích přitažlivosti dvou závaží vznikají při obarvení fraktální vzory (ve kterých se kupodivu nachází i oblasti přitažlivosti nejvzdálenějšího magnetu).

fractals22_3
Obrázek 3: Detailní pohled na dělení komplexní plochy do třech oblastí

4. Odvození funkce Newtonovy iterační metody pro polynom z3-1=0

Při odvozování funkce použité v Newtonově iterační metodě budeme vycházet z výše uvedeného vztahu:

zi+1=zi-f(zi)/f'(zi)

Do tohoto vztahu dosadíme polynom z3-1, tj. bude platit:

f(z)=z3-1
f'(z)=3z2

A po dosazení do obecného iteračního vztahu:

zi+1=zi-(zi3-1)/(3zi2)=
=zi(3zi2)/(3z­i2)-(zi3-1)/(3zi2)=
=2zi3/(3zi2)+­1/(3zi2)=
=2/3zi+1/(3z2

Po úpravách a následném rozpisu posledního výrazu na jeho reálnou a imaginární složku dostaneme výsledek vhodný pro programové zpracování:

Re: zre i+1 = (2/3)zre i + (zre i2-zim i2)/(3(zre i4+zim i4+2zre i2 zim i2))
Im: zim i+1 = (2/3)zim i – (2zre i zim i)/(3(zre i4+zim i4+2zre i2 zim i2))

5. Newtonův fraktál polynomu z3-1

Výše uvedené výrazy pro reálnou a imaginární složku hodnoty zi+1 jsou použity ve funkci recalcNewton(), která provádí výpočet a následné zobrazení Newtonovy fraktální množiny. Tato funkce pracuje následujícím způsobem: nejprve je každý bod v komplexní rovině chápán jako počáteční hodnota z0, tj. hodnota, pro kterou má být nalezen kořen polynomu. Následně je tato hodnota postupně iterována, přičemž iterační smyčka skončí v případě, že se hodnota zk (pro k>=0) dostatečně přiblížila k jednomu ze tří kořenů polynomu zi+1. Bod je poté obarven jednou ze tří barev podle toho, ke kterému kořenu řešení dospělo. Iterační smyčka se však také ukončí v případě, že bylo dosaženo maximálního množství iterací bez toho, aby se řešení k nějakému kořenu dostalo – v tomto případě je bodu přiřazena černá barva. Toto alternativní ukončení iterační smyčky je nutné v programu ponechat, neboť pro některé počáteční hodnoty z0 může výpočet trvat i nekonečně dlouhou dobu (například se jedná o počátek, tj. bod 0+0i). Kód popsané funkce může v programovacím jazyku C vypadat následovně:

//-----------------------------------------------------------------------------
// Překreslení Newtonovy fraktální množiny pro polynom z^3-1=0
//-----------------------------------------------------------------------------
void recalcNewton(pixmap *pix,                  // pixmapa pro vykreslování
                  int    maxiter,               // maximální počet iterací
                  double scale,                 // měřítko obrazce
                  double xpos,                  // posun obrazce
                  double ypos)
{
    double zx, zy, zx2, zy2, zxn, zyn;          // složky komplexní proměnné Z a Z^2
    double zx0, zy0;
    double xmin, ymin, xmax, ymax;              // rohy vykreslovaného obrazce
                                                // v komplexní rovině
    int    x, y;                                // počitadla sloupců a řádků v pixmapě
    int    iter;                                // počitadlo iterací
    unsigned char r, g, b;

#define PI 3.1415927
#define EPSILON 0.1
    // makro vyjadřující druhou mocninu vzdáleností dvou bodů
#define DIST2(x1, y1, x2, y2) (((x1)-(x2))*((x1)-(x2))+((y1)-(y2))*((y1)-(y2)))

    // kořeny polynomu z^3-1=0 můžeme dopředu přesně vypočítat
    double rootx[3]={1.0, cos(2.0*PI/3.0), cos(4.0*PI/3.0)};
    double rooty[3]={0.0, sin(2.0*PI/3.0), sin(4.0*PI/3.0)};

    calcCorner(xpos, ypos, scale, &xmin, &ymin, &xmax, &ymax);
    zy0=ymin;
    for (y=0; y<pix->height; y++) {             // pro všechny řádky v pixmapě
        zx0=xmin;
        for (x=0; x<pix->width; x++) {          // pro všechny pixely na řádku
            zx=zx0;                             // nastavit počáteční hodnotu Z(0)
            zy=zy0;
            for (iter=0; iter<maxiter; iter++) {// iterační smyčka
                zx2=zx*zx;                      // zkrácený výpočet druhé mocniny složek Z
                zy2=zy*zy;
                zxn=2.0/3.0*zx+(zx2-zy2)/(3.0*(zx2*zx2+zy2*zy2+2.0*zx2*zy2));
                zyn=2.0/3.0*zy-2.0*zx*zy/(3.0*(zx2*zx2+zy2*zy2+2.0*zx2*zy2));
                                                // test na přiblížení ke kořenům polynomu
                if (DIST2(zxn, zyn, rootx[0], rooty[0])<EPSILON ||
                    DIST2(zxn, zyn, rootx[1], rooty[1])<EPSILON ||
                    DIST2(zxn, zyn, rootx[2], rooty[2])<EPSILON) break;
                zx=zxn;
                zy=zyn;
            }
            r=g=b=0x00;                         // černá barva je nastavena v případě
                                                // že není rozhodnuto, ke kterému
                                                // kořenu bod přísluší
            if (DIST2(zxn, zyn, rootx[0], rooty[0])<EPSILON) r=0xff, g=0x00, b=0x00;
            if (DIST2(zxn, zyn, rootx[1], rooty[1])<EPSILON) r=0x00, g=0xff, b=0x00;
            if (DIST2(zxn, zyn, rootx[2], rooty[2])<EPSILON) r=0x00, g=0x00, b=0xff;
            putpixel(pix, x, y, r, g, b);
            zx0+=(xmax-xmin)/pix->width;        // posun na další bod na řádku
        }
        zy0+=(ymax-ymin)/pix->height;           // posun na další řádek
    }
} 

Funkce recalcNewton() je použita v dnešním demonstračním příkladu, ve kterém je možné po překladu a následném spuštění zobrazit Newtonovu fraktální množinu pro polynom z3-1=0. Na obrázcích 4–7 je zobrazena vzniklá Newtonova množina pro různý počet iterací. Všimněte si, že pokud je provedena pouze jedna iterace, přísluší ke třem kořenům pouze body ze tří velkých oblastí (kořeny leží přibližně ve středech těchto oblastí) a dále pak tři mnohem menší oblasti, které vůči svým větším protějškům leží na opačné straně počátku – to vychází z principu výpočtu, kdy se komplexní číslo při umocnění otáčí okolo počátku.

fractals22_4
Obrázek 4: Newtonova fraktální množina po jedné iteraci

Pokud jsou při výpočtu Newtonovy fraktální množiny použity maximálně dvě iterace, je v celé komplexní rovině vytvořeno už patnáct obarvených oblastí, přičemž je opět patrná symetrie vůči počátku – to je patrné také z pátého obrázku. S rostoucím počtem iterací (viz obrázky šest a sedm) se komplexní rovina stále více fragmentuje na barevné oblasti, na druhou stranu se tři největší oblasti (ve kterých je řešení polynomu stabilní) neustále zvětšují a rostou na Riemannově sféře směrem k nekonečnu. Pokud bychom provedli nekonečné množství iterací, byla by prakticky celá komplexní rovina rozdělena obarvena třemi barvami, části obarvené černě (tj. „nerozhodný stav“) by měly nulovou míru.

fractals22_5
Obrázek 5: Newtonova fraktální množina po dvou iteracích

fractals22_6
Obrázek 6: Newtonova fraktální množina po deseti iteracích

fractals22_7
Obrázek 7: Newtonova fraktální množina po 255 iteracích

6. Modifikovaný Newtonův fraktál polynomu z3-1=0

Na předchozích čtyřech obrázcích byl Newtonův fraktál zobrazen tím způsobem, že jednotlivé body byly obarveny na základě kořenu, ke kterému iterační výpočet směřoval. Tento fraktál je však možné zobrazit i modifikovaným způsobem, který vychází z podobné myšlenky, jakou jsme použili už při zobrazování Mandelbrotových a Juliových množin. Jednotlivé body nejsou obarveny na základě příslušnosti k některému z kořenů, ale podle počtu iterací, které bylo nutné provést pro jednoznačné rozhodnutí, ke kterému kořenu bod přísluší. Proceduru, která provádí výpočet takto modifikovaného Newtonova fraktálu, si ukážeme v dalším pokračování tohoto seriálu; detail tímto způsobem vizualizované Newtonovy množiny je zobrazen na osmém obrázku, ve kterém je počet iterací vyjádřen stupněm šedé barvy, kde černá barva odpovídá pouze jedné iteraci a bílá barva maximálnímu počtu iterací.

fractals22_8
Obrázek 8: Body v Newtonově fraktální množině zobrazené na základě počtu iterací

7. Newtonovy fraktály polynomů vyšších stupňů

Již z popisu tvorby Newtonovy fraktální množiny na základě polynomu z3-1=0 je zřejmé, že se nemusíme omezovat pouze na třetí mocninu hodnoty z. Skutečně, je možné použít buď libovolnou celočíselnou mocninu větší než 2 (potom je výsledný obrázek symetrický), libovolnou reálnou mocninu (symetrie je narušena) a dokonce se může vztah upravit tak, že se za mocninu dosadí komplexní číslo – v tom případě je výsledný obraz zkreslen ještě více, než v případě reálné mocniny. Bližší informace o těchto modifikacích Newtonovy fraktální množiny si řekneme v dalším pokračování tohoto seriálu.

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

V následujících dvou částech tohoto seriálu si ukážeme tvorbu modifikované Newtonovy množiny a také si popíšeme dva nové typy fraktálů vytvořených v komplexní rovině. Jedná se o fraktál pojmenovaný Magnet a o fraktál pojmenovaný Barnsley.

Našli jste v článku chybu?

21. 3. 2006 16:08

My jsme ruzne iteracni metody pro hledani korenu na vysce taky probirali. Bylo to vsak v numericke matematice, kde se za jeden semestr musela probrat spousta dalsich veci, takze jsme o ni moc nediskutovali. Je to skoda, stejne tak jsme nestihli i dalsi zajimava temata, napriklad aproximacni a interpolacni krivky (ty samozrejme znam z grafiky, matematici vsak pouzivaji i dalsi typy krivek, zejmena jsou zajimave ty odvozovacky, ktere jsou zalozeny na trosku jinem principu nez ty "graficke&quo…

21. 3. 2013 7:39

Jožo (neregistrovaný)

...fraktály sů super... :)

Měšec.cz: Finančním poradcům hrozí vracení provizí

Finančním poradcům hrozí vracení provizí

Vitalia.cz: Říká amoleta - a myslí palačinka

Říká amoleta - a myslí palačinka

Podnikatel.cz: Přehledná titulka, průvodci, responzivita

Přehledná titulka, průvodci, responzivita

120na80.cz: Rakovina oka. Jak ji poznáte?

Rakovina oka. Jak ji poznáte?

DigiZone.cz: Česká televize mění schéma ČT :D

Česká televize mění schéma ČT :D

Lupa.cz: Avast po spojení s AVG propustí 700 lidí

Avast po spojení s AVG propustí 700 lidí

DigiZone.cz: Sony KD-55XD8005 s Android 6.0

Sony KD-55XD8005 s Android 6.0

Vitalia.cz: Láska na vozíku: Přitažliví jsme pro tzv. pečovatelky

Láska na vozíku: Přitažliví jsme pro tzv. pečovatelky

Lupa.cz: Proč firmy málo chrání data? Chovají se logicky

Proč firmy málo chrání data? Chovají se logicky

Vitalia.cz: Jsou čajové sáčky toxické?

Jsou čajové sáčky toxické?

Lupa.cz: Babiš: E-shopů se EET možná nebude týkat

Babiš: E-shopů se EET možná nebude týkat

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

Recenze Westworld: zavraždit a...

Podnikatel.cz: EET: Totálně nezvládli metodologii projektu

EET: Totálně nezvládli metodologii projektu

Vitalia.cz: Chtějí si léčit kvasinky. Lék je jen v Německu

Chtějí si léčit kvasinky. Lék je jen v Německu

Vitalia.cz: To není kašel! Správná diagnóza zachrání život

To není kašel! Správná diagnóza zachrání život

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

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

Lupa.cz: UX přestává pro firmy být magie

UX přestává pro firmy být magie

Vitalia.cz: Baletky propagují zdravotní superpostel

Baletky propagují zdravotní superpostel

Podnikatel.cz: Víme první výsledky doby odezvy #EET

Víme první výsledky doby odezvy #EET

120na80.cz: Jak oddálit Alzheimera?

Jak oddálit Alzheimera?