2. Domácí úkol: Newtonův fraktál

Úvod

Cílem druhého úkolu je dále prozkoumat vlastní typy, práci s maticemi a poli a v neposlední řadě i s vizualizačními nástroji dostupnými v Julia.

Obrázek 1. Newtonův fraktál pro funkci $f(z) = z^3 - 1$ , která má tři komplexní kořeny.

Teoretický background

Newtonova metoda, také často označovaná jako Newtonova-Raphsonova metoda, má komplikovanou historii. V jisté formě ji zřejmě dokázali používat už Babyloňané, její moderní formulace pravděpodobně pochází až ze 17. století.

Dnes ji můžeme popsat snadno takto: pro komplexní funkci komplexní proměnné $f: \mathbb{C} \to \mathbb{C}$ a zadané komplexní číslo $z_0 \in \mathbb{C}$ (počáteční aproximace, počáteční podmínka) sestrojme rekurentní posloupnost

z_{n+1} = z_n - \frac{f(z_n)}{f'(z_n)}, \quad n = 0,1,2,\ldots,

kde $f'$ označuje derivaci funkce $f$ . Naším přáním je, aby takto zkonstruovaná posloupnost konvergovala k jistému bodu $z_*$ , který je nulovým bodem funkce $f$ , tj. $f(z_*) = 0$ .

Pokud je první aproximace (tj. $z_0$ ) blízko u nulového bodu funkce $f$ , pak lze takovéto chování očekávat. Ukazuje se, že v případě bodů vzdálenějších od nulových bodů je chování rekurentní posloupnosti popsané výše daleko zajímavější. Komplexní rovina se nerozdělí na souvislé oblasti konvergence, ale takovéto rozdělení vykazuje fraktální kvality. Ve skutečnosti není nijak zaručeno, že uvedená rekurentní posloupnost bude vždy konvergentní, nebo že ji vždy lze zkonstruovat (ve jmenovateli zlomku se může vyskytnout nula, můžeme se dostat mimo definiční obor funkce $f$ , atp.).

Naším cílem v tomto úkolu je tuto metodu implementovat a graficky prozkoumat její chování.

Implementační pokyny

Své řešení opět doplňte do souboru main.jl, který obsahuje hrubou kostru řešení. Odpovídající testy naleznete v adresáři tests (můžete přidat případně i své vlastní).

Typ `Newton{T}`

Objekty tohoto typu popisují obdélníkovou mřížku v komplexní rovině určenou rozsahy re_min a re_max v reálné části a im_min a im_max v imaginární části. V reálném směru máme na mřížce width bodů a podobně v imaginárním směru height bodů. Komplexní čísla reprezentujeme typem Complex{T}, kde T je podtyp typu AbstractFloat.

HINT: Kdyby například re_min = -1.0, re_max = 1.0, im_min = -1.0 a im_max = 1.0 a dále width = height = 3, pak uvažujeme body -1.0-1.0im, -1.0+0.0im, -1.0+1.0im, 0.0-1.0im, 0.0+0.0im, 0.0+1.0im, 1.0-1.0im, 1.0+0.0im a 1.0+1.0im.

Dále objekty tohoto typu nesou informaci o sledované funkci f a její derivace df. Tyto funkce berou jeden (komplexní) argument a vrací jednu (komplexní) hodnotu. Je na uživateli, aby zadal správné funkce, jejich vztah (jedna je derivací druhé) nijak testovat nebudeme.

V neposlední řadě tyto objekty obsahují informaci o napočtených experimentech. Počáteční podmínky pro iterování volíme jako body na uvedené mřížce a získané informace budeme zaznamenávat do matic (viz dále, které indexujeme obvyklým způsobem: bodu vlevo dole (tj. re_min + im_min * 1im) odpovídá index [1, 1] a bodu vpravo nahoře (tj. re_max + im_max * 1im) odpovídá index [width, height]. Zaznamenáváme následující:

matice data obsahuje limity (k detekci konvergence viz dále) posloupností s odpovídajícími počátečními podmínkami, nebo komplexní NaN v případě nedosažení konvergence,
matice iterations obsahuje provedený počet (Int64) iterací pro odpovídající počáteční podmínku, případně 0 pokud během iterování došlo k přetečení nebo se během iterování vyskytl komplexní NaN, nebo nastala výjimka.
matice roots obsahuje příslušnost odpovídající počáteční podmínky k hypotetickému nulovému bodu zkoumané funkce (ty indexujeme od 1), případně 0 pokud pro tuto počáteční podmínku nedošlo k dosažení konvergence.

Konečně vektor roots_list komplexní čísel obsahuje přibližné hodnoty nalezených nulových bodů (kořenů) uvažované funkce.

Pokud to bude potřeba, můžete si přidat i další pomocné či pracovní atributy.

Vnitřní konstruktor tohoto typu očekává popořadě funkci, její derivaci, rozsahy reálných a komplexních složek. Rozměry mřižky lze předat jako keyword argumenty a výchozí hodnota je nastavena na $500 \times 500$ mřížku. Tento konstruktor připraví místo pro výpočty (matice vhodných rozměrů, atd.), ale ještě automaticky nespustí výpočet. Dále ověří, že rozsahy dávají smysl (např. re_min není větší nebo rovno re_max) a že rozměry mřížky v každém směru jsou alespoň $2$ . V případě nesplnění těchto podmínek vyvolá ArgumentError.

Dále je vhodné pro tento nový typ implementovat show metodu, jinak bude výpis objektu nepřehledný, zde nechávám volnost vaší kreativitě.

Metoda `newton`

Tato metoda je ústřední pro naše snažení a implementuje samotnou Newtonovu metodu. Přijímá popořade poziční argumenty

f a df představující funkci a její derivaci, ty jsou schopny pracovat s komplexní proměnnou a vrací komplexní číslo,
komplexní číslo z představující počáteční podmínku (komplexní číslo typu S <: Complex),

a dále keyword argumenty

maxiter představující maximální povolený počet iterací (výchozí hodnota 1_000),
atol je přesnost použitá pro rozhodnutí o konvergenci (viz dále; výchozí hodnota je sqrt(eps(abs(one(S))))).

Metoda implementuje Newtonovu metodu popsanou výše a ukončuje se následovně:

Pokud se od sebe dvě poslední iterace v absolutní hodnotě neliší více než o atol, pak končíme úspěchem a vracíme dvojici, kde v první složce je hodnota poslední iterace a v druhé provedený počet iterací.
Pokud dosáhneme maximálního počtu iterací, pak vracíme dvojici, kde v první složce je S(NaN) a v druhé onen maximální počet iterací.
Pokud během výpočtu dojde v rámci aritmetických operací ve strojových číslech k vytvoření NaN nebo Inf hodnoty, nebo dojde k problému s vyhodnocováním funkčních hodnot zadaných funkcí, pak výpočet zcela selhal a vracíme dvojici, kde v první složce je S(NaN) a v druhé 0.

HINT: V kódu je možné pro testování různých NaN hodnot použít metodu isnan, podobně existuje isinf metoda.

Metoda `compute!`

Tato metoda přijímá objekt nt typu Newton{T} a jako keyword argumenty opět maximální počet iterací maxiter (výchozí hodnota 1_000) a absolutní toleranci atol. Tato metoda provádí finální výpočet ve dvou krocích:

V každém bodě mřížky provede Newtonovu metodu a výsledek zapíše do odpovídajícího prvku v maticích data a iterations (můžete přirozeně přímo použít metodu newton zmíněnou výše, ale není to nutné).
Dále odhalí nalezené různé limitní body a to následujícím elementárním způsobem: nalezené kořeny jsou přidávány do pole roots_list, při procházení matice data ignorujeme NaN hodnoty, pokud se daný prvek matice data liší od některého z již nalezených kořenů o méně než sqrt(atol), pak je považujeme za totožné a do odpovídajícího prvku matice roots vložíme index tohoto kořene (jakožto prvku pole roots_list). Pokud uvažovaný prvek matice data není dostatečně blízko žádnému z prvků pole roots_list, tak ho vložíme na konec tohoto pole jakožto nový kořen a do odpovídajícího prvku matice roots opět zaznamenáme informaci o indexu.

Metoda vrací zpracovávaný objekt nt.

Metoda `show`

Zde vám ponechávám největší volnost. Metoda show přijme objekt typu Newton{T} a vytvoří jeho grafickou reprezentaci pomocí vámi zvoleného grafického balíčku (Images, PyPlot, atd., viz materiály). Pro ilustraci přikládám obrázek na začátku tohoto zadání.

Typicky se při vizualizaci pro každý nalezený kořen volí odstín barvy a počet iterací pak modifikuje sytost barvy. K tomuto účelu se velmi hodí způsob zadání barvy pomocí HSV (hue, saturation, value). Ale nebojte se zkusit data obarvit dle své fantazie.

Pozor, tento krok pravděpodobně naruší testovací prostředí. Je potřeba do něj příslušný balíček doinstalovat editací .gitlab-ci.yml souboru, kam před spuštěním testů přidáme řádek tvaru

    - julia -e 'import Pkg; Pkg.add("Images")'

kde místo Images uvedete vámi použitý balíček. Pokud s tímto krokem budete bojovat, tak stačí odevzdat i řešení, kde tento balíček nebudete explicitně používat (zakomentování using …) a při kontrole si s tím poradím lokálně. Uveďte tuto informaci v textu svého MR.

Tuto metodu vzhledem k její otevřenosti žádným způsobem netestuji v přiložených testech.

Ukázka použití této implementace je v notebooku examples.ipynb. Zkuste se jím inspirovat a vytvořit vlastní zajímavý obrázek (jiná funkce, jiné barevné znázornění, atp.). Z nejhezčích výtvorů vytvořím plakát pro zkrášlení šedých fakultních chodeb.

Spouštění testů

Zadání je doplněno testy cílícími na vlastnosti typu Newton{T} a metod newton a compute!.

$ julia test/runtests.jl
┌ Info:
│ Running tests in test_newton.jl...
└ ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
Test Summary:                                                | Pass  Total  Time
Newton: constructor                                          |   13     13  0.2s
Test Summary:                                                | Pass  Total  Time
newton: bad cases                                            |   10     10  0.1s
Test Summary:                                                | Pass  Total  Time
newton: unique roots                                         |    3      3  0.0s
Test Summary:                                                | Pass  Total  Time
newton: multiple roots, starting nearby                      |    4      4  0.1s
Test Summary:                                                | Pass  Total  Time
newton: iteration limit                                      |    2      2  0.0s
Test Summary:                                                | Pass  Total  Time
compute!: one root case                                      |    7      7  0.4s
Test Summary:                                                | Pass  Total  Time
compute!: two roots                                          |  103    103  0.4s
Test Summary:                                                | Pass  Total  Time
compute!: division by zero                                   |    3      3  0.0s

Tyto testy se spouští i automaticky na Gitlabu a na stránce Merge Requestu jejich výsledek uvidíte ve formě červené nebo zelené ikonky.

Součástí vyhodnocení řešení bude porovnání rychlosti i robustnosti vaší implementace vzhledem k ostatními řešením.

Varování:

Pozor, pokud přidáváte grafických balíček pomocí using, tak ho přidejte do .gitlab-ci.yml, viz poznámku výše u metodu show.

Řešení a odevzdání

Opět vytvořte větev odvozenou z větve assignment/02-newton a nezvěte ji solution/02-newton. Do solution/02-newton vložte své řešení editací souboru main.jl případně přidáním testů do složky test. Až budete se svým řešením spokojeni, vytvořte MR (to můžete i dříve, aspoň uvidíte výsledek testů, pokud je nespouštíte lokálně) s cílovou větví assignment/02-newton a přiřaďte mě k němu jako assignee. Tímto aktem úkol odevzdáte.

2. Domácí úkol: Newtonův fraktál

Úvod

Teoretický background

Implementační pokyny

Typ Newton{T}

Metoda newton

Metoda compute!

Metoda show