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

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 z³-1=0
5. Newtonův fraktál polynomu z³-1=0
6. Modifikovaný Newtonův fraktál polynomu z³-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 a_nxⁿ+a_n-1x^n-1…a₂x²+a₁x+a₀ 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:

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

přičemž je levá strana rovnice (ve výše uvedeném vzorci označena jako x_i+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ě x₀, 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ě.

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 x²-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 x²+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 i²=-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ě:

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

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ě z₀. 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 z³-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 z₀ 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 z₀ 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.

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 z_0A a z_0B 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).

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

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

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

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

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

f(z)=z³-1
f'(z)=3z²

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

z_i+1=z_i-(z_i³-1)/(3z_i²)=
=z_i(3z_i²)/(3z­_i²)-(z_i³-1)/(3z_i²)=
=2z_i³/(3z_i²)+­1/(3z_i²)=
=2/3z_i+1/(3z_i­²) 

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: z_{re i+1} = (2/3)z_{re i} + (z_{re i}²-z_{im i}²)/(3(z_{re i}⁴+z_{im i}⁴+2z_{re i}² z_{im i}²))
Im: z_{im i+1} = (2/3)z_{im i} – (2z_{re i} z_{im i})/(3(z_{re i}⁴+z_{im i}⁴+2z_{re i}² z_{im i}²))

5. Newtonův fraktál polynomu z³-1

Výše uvedené výrazy pro reálnou a imaginární složku hodnoty z_i+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 z₀, 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 z_k (pro k>=0) dostatečně přiblížila k jednomu ze tří kořenů polynomu z_i+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 z₀ 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 z³-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.

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.

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

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

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

6. Modifikovaný Newtonův fraktál polynomu z³-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í.

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 z³-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.