Newtoni meetod ja fraktaalsed basseiniääred

Varasemates näidetes oleme näinud, et diskreetse või pideva dünaamika atraktorid võivad olla fraktaalsed. Antud tööleht näitab, et atraktiivsuse basseinide ääred samuti võivad olla fraktaalsed. Me uurime neid konkreetse näite puhul, mis on Newtoni meetod funktsioonide nullkohtade leidmiseks. Siin kasutame omadust, et nullkohad on atraktorid, ja uurime nende atraktiivsuse basseine. Nende uurimiseks kasutame SymPy ja NumPy.

In [1]:
from PIL import Image
import IPython.display as disp
import numpy as np
import sympy as sp
from io import BytesIO
import multiprocessing as mp

Kuna joonistused on suuremad kui tavaliselt, suurendame ka tulemusala.

In [2]:
%%html
<style>
.output_wrapper, .output {
    height:auto !important;
    max-height:10000px;
}
.output_scroll {
    box-shadow:none !important;
    webkit-box-shadow:none !important;
}
</style>

Atraktorid ja nende atraktiivsuse basseinid

Tuletame meelde atraktori definitsiooni:

Faasiruumi alamhulka $A \subset Q$ nimetatakse atraktoriks, kui kehtib:

  • Igas punktis $x \in A$ kehtib $f(x) \in A$.
  • On olemas lahtine ümbrus $\mathcal{B}(A) \supset A$, mida nimetatakse $A$ atraktiivsuse basseiniks, kust dünaamika koondub atraktori juurde. See tähendab et $\mathcal{B}(A)$ sisaldab neid punkte $x \in Q$ omadusega, et iga lahtise ümbruse $N \supset A$ korral leidub $n \in \mathbb{N}$ nii, et iga $m \geq n$ korral $f^m(x) \in N$. Atraktiivsuse bassein on ehk hulk $$\mathcal{B}(A) = \{x \in Q \,|\, \forall N \supset A: \exists n \in \mathbb{N}: \forall m \geq n: f^m(x) \in N\}.$$
  • Ei ole võimalik hulka $A$ jagada mitmesse ossa, mis rahuldavad neid tingimusi.

Me pöörame nüüd tehelepanu atraktiivsuse basseinile. Kuna $\mathcal{B}(A)$ on lahtine hulk, järeldub et tema Hausdorff-Besicovitši ja Minkowski-Bouligandi dimensioonid on samad kui faasiruumi $Q$ dimensioon $\dim Q$. Naiivselt võib arvata, et ääre $\partial\mathcal{B}(A)$ dimensioon on $\dim Q - 1$, aga selline järeldus ei kehti.

Newtoni meetod

Varasemas loengus kasutasime Newtoni meetodi võrrandite numbriliseks lahendamiseks. Nüüd uurime Newtoni meetodu kui dünaamilist süsteemi. Selleks eeldame, et meie faasiruum oleks kompleksarvude hulk $Q = \mathbb{C}$, ja meil on antud kaks korda diferentseeruv funktsioon $f: \mathbb{C} \to \mathbb{C}$, millel on vähemalt üks nullkoht $z_*$, kus kehtib $f(z_*) = 0$ ja $f'(z_*) \neq 0$. Nullkoha leidmiseks valitakse algtingimust $z_0 \in \mathbb{C}$ ja arvutatakse trajektoori dünaamika

$$z_{n + 1} = F(z_n) = z_n - \frac{f(z_n)}{f'(z_n)}$$

alusel. Funktsiooni $f$ nullkoht $z_*$ on süsteemi püsipunkt, kuna kehtib

$$F(z_*) = z_* - \frac{f(z_*)}{f'(z_*)} = z_*.$$

Selle püsipunkti stabiilsuse määramiseks arvutame dünaamika tuletist

$$F'(z) = 1 - \frac{f'(z)^2 - f(z)f''(z)}{f'(z)^2} = \frac{f(z)f''(z)}{f'(z)^2}.$$

Nullkohas $z_*$ võrdub see nulliga. Sellest järeldub, et $z_*$ on atraktor, ja seega on olemas ka atraktiivsuse bassein.

Järgmisena uurime, millisesse basseini punktid $z \in \mathbb{C}$ kuuluvad. Selleks värvime iga basseini unikaalse värviga, mida me defineerime järgnevalt.

In [3]:
def color(n, i):
    k = int((1530 * i) / n)
    return (max(0, min(255, 510 - abs(k - 1020))), max(0, min(255, 510 - abs(k - 510))), max(0, min(255, abs(k - 765) - 255)))

Newtoni meetodi paralleelseks rakendamiseks defineerime klassi, mida me saame mitme protsessori peal korraga käivitada. Klassi konstruktori parameetrid on järgnevalt:

  • expr on funktsioon $f(z)$, SymPy string formaadis.
  • zero on funktsiooni nullkohad $z_{*1}, \ldots, z_{*p}$, mida me soovime värvidega tähistada.
  • eps on kaugus $\epsilon$: nullkoht on leitud, kui $|z_n - z_{*k}| < \epsilon$.
  • maxit on iteratsioonisammude ülempiir.

Klassiobjekti võib käivitada ühe parameetriga z0, mis on algtingimus $z_0$.

In [4]:
class newtit(object):
    def __init__(self, expr, zero, eps, maxit):
        # Jätame parameetrid meelde.
        self.expr = expr
        self.zero = zero
        self.n0 = len(zero)
        self.eps = eps
        self.maxit = maxit
        self.f = None
        self.df = None
        self.nf = None
        
    def __call__(self, z0):
        # Esimesel käivitamisel arvutame f'(z), konverteerime numbriliseks arvutamiseks ja jätame meelde.
        if self.f == None:
            z = sp.Symbol('z')
            fkt = sp.sympify(self.expr)
            dfkt = sp.diff(fkt, z)
            nwt = sp.simplify(z - fkt / dfkt)

            self.f = sp.lambdify(z, fkt, 'numpy')
            self.df = sp.lambdify(z, dfkt, 'numpy')
            self.nf = sp.lambdify(z, nwt, 'numpy')

        zn = z0
        for n in range(self.maxit):
            # Kui f'(z) = 0, tähistame punkti musta värviga.
            if self.df(zn) == 0:
                return (0, 0, 0)
            # Newtoni samm.
            zn = self.nf(zn)
            # Kas me oleme ühe nullkoha lähedal?
            dist = [abs(z - zn) for z in self.zero]
            dm = min(dist)
            if dm < self.eps:
                k = dist.index(dm)
                return color(self.n0, k)
        # Kui trajektoor ei koondu, tähistame punkti musta värviga.
        return (0, 0, 0)

Lisaks defineerime funktsiooni, mis rakendab Newtoni meetodi punktidele ristkülikukujulises piirkonnas ja joonistab tulemusi.

In [5]:
def newton(expr, zero, zmin, zmax, nx, ny, eps, maxit):
    # Ristküliku piirid.
    xmin = zmin.real
    ymin = zmin.imag
    xmax = zmax.real
    ymax = zmax.imag

    # Alguspunktid z0, millele me soovime Newtoni meetodi rakendada.
    z0 = [complex(xmin + i * (xmax - xmin) / (nx - 1), ymax - j * (ymax - ymin) / (ny - 1)) for j in range(ny) for i in range(nx)]
    
    # Rakendame funktsiooni paralleelselt.
    with mp.Pool() as pool:
        res = pool.map(newtit(expr, zero, eps, maxit), z0)
    
    # Joonistamiseks konverteerime tulemusi graafikaformaati.
    col = np.array(res, dtype=np.uint8).reshape(ny, nx, 3)
    b = BytesIO()
    imr = Image.fromarray(col)
    imr.save(b, format='png')
    
    # Näitame tulemust.
    disp.display_png(b.getvalue(), raw=True)

Eelnevalt defineeritud funktsioonid saame nüüd rakendada konkreetsetele näidetele.

Näide: $z^p - 1 = 0$

Atraktiivsuse basseinide uurimiseks vaatame näidet $f(z) = z^p - 1$, millel on $p$ nullkohta ühikringi $|z| = 1$ peal ja rahuldavad

$$z_{*k} = e^{\frac{2\pi ik}{p}}.$$

Selleks defineerime funktsiooni:

In [6]:
def newtpow(p, zmin, zmax, nx, ny, eps, maxit):
    expr = "z**%d - 1" % p
    zero = [np.exp(k * 2.0j * np.pi / p) for k in range(p)]
    newton(expr, zero, zmin, zmax, nx, ny, eps, maxit)

Nüüd joonistame tulemusi erinevate parameetri $n$ väärtuste puhul:

In [7]:
for p in range(3, 11):
    newtpow(p, -2-2j, 2+2j, 1024, 1024, 1e-3 / p, 1024)

Jonistused näitavad, et igal nullkohal $z_*$ on atraktiivsuse bassein, ja et basseini äärtel on fraktaalne struktuur. Selle struktuuri täpsemaks uurimiseks suurendame osa $[-0.8, 0] \times [-0.3i, 0.3i]$ pildist, mis vastab parameetri väärtusele $p = 3$:

In [8]:
newtpow(3, -0.8-0.3j, 0.3j, 1024, 768, 1e-6, 1024)

Parema serva keskel on punkt $z = 0$, mis ei kuulu ühesse basseini, vaid basseinide äärde, kuna ta asub just basseinide puutepunktis. See tähendab, et Newtoni meetod selles punktis ei koondu. Sellest on lihtne aru saada, kuna seal kehtib $f'(z) = 3z^2 = 0$, ja vastav trajektoor hajub. Kuna tuletisel $f'(z)$ on ainult üks nullkoht $z = 0$, ehk hulgal

$$C_0 = \{z \in \mathbb{C}, f'(z) = 0\} = \{0\}$$

on ainult üks liige, leiame, et trajektoor $(z_n, n \in \mathbb{N})$ hajub siis kui on olemas $n \in \mathbb{N}$ nii, et $z_n = 0$. Trajektoori leiame kujutuse

$$F(z) = z - \frac{f(z)}{f'(z)} = \frac{1}{3z^2} + \frac{2z}{3}$$

abil. Algtingimused, millele sellised trajektoorid vastavad, on näiteks sellised:

  • $n = 1$: Me leiame, et $z_1 = 0$ siis kui $$z_0 \in C_1 = \{z \in \mathbb{C}, F(z) \in C_0\} = \left\{-\frac{1}{\sqrt[3]{2}}, \frac{e^{\frac{\pi}{3}i}}{\sqrt[3]{2}}, \frac{e^{-\frac{\pi}{3}i}}{\sqrt[3]{2}}\right\}.$$ Esimene punkt $z_0 = -\frac{1}{\sqrt[3]{2}}$ on joonisel vasaku serva lähedal, kus basseinide tipud kohtuvad.
  • $n = 2$: On lihtne näha et $z_2 = 0$ siis kui $z_1 \in C_1$. Samuti on olemas hulk $$C_2 = \left\{z \in \mathbb{C}, F(z) \in C_1\right\},$$ milles on 9 punkti.
  • Sama printsiibiga võib defineerida hulke $C_3, C_4, \ldots$ üldise valemi $$C_{n + 1} = \left\{z \in \mathbb{C}, F(z) \in C_n\right\}$$ järgi, mis on samuti basseini ääre alamhulgad. Nende mahukused on $\#C_n = 3^n$, mis näitab et basseinte ääre mahukus on vähemalt loenduv.

Lõpuks tasub mainida, et basseinte ääre mahukus on tegelikult mitteloenduv - nagu ka Cantori hulk ei ole loenduv. Seda saab tõestada sama meetodiga. Fraktaali omadust saab näidata näiteks Minkowski-Bouligandi dimensiooni arvutamisel. Sama printsiip kehtib ka teiste $p$ väärtuste puhul, näiteks $p = 4$:

In [9]:
newtpow(4, 0, 0.75+0.75j, 1024, 1024, 1e-6, 1024)

Tuletisest $f'(z) = 4z^3$ järeldub, et trajektoor hajub, kui ta sisaldab punkti $z = 0$. See punkt asub pildil vasakul alumises nurgas. Seega defineerime ka sel juhul:

$$C_0 = \{z \in \mathbb{C}, f'(z) = 0\} = \{0\}.$$

Dünaamilise süsteemi kujutuse

$$F(z) = z - \frac{f(z)}{f'(z)} = \frac{1}{4z^3} + \frac{3z}{4}$$

abil leiame veel punkte basseini ääres. Esimene samm selleks on hulk

$$C_1 = \{z \in \mathbb{C}, F(z) \in C_0\} = \left\{\frac{1 + i}{\sqrt{2}\sqrt[4]{3}}, \frac{1 - i}{\sqrt{2}\sqrt[4]{3}}, \frac{-1 + i}{\sqrt{2}\sqrt[4]{3}}, \frac{-1 - i}{\sqrt{2}\sqrt[4]{3}}\right\}.$$

Pildil asub esimene punkt $z = \frac{1 + i}{\sqrt{2}\sqrt[4]{3}}$ "tärni" keskel, basseinide puutepunktis.

Näide: $\prod_{k = 1}^p(z - z_{*k})$

Me võime oma tulemusi üldistada, kui me valime suvalise polünoomi nullkohtadega $z_1, \ldots, z_p$, see tähendab et

$$f(z) = \prod_{k = 1}^p(z - z_{*k}).$$

Lihtsamaks joonistamiseks defineerime sellist funktsiooni:

In [10]:
def newtpoly(zero, zmin, zmax, nx, ny, eps, maxit):
    expr = "*".join(["(z - (%s))" % z for z in zero])
    newton(expr, zero, zmin, zmax, nx, ny, eps, maxit)

Kui me valime näiteks

$$f(z) = z(z - 1)(z + 1)(z - i) (z + i) = z^5 - z,$$

siis tekib selline pilt:

In [11]:
newtpoly([0, 1, -1, 1j, -1j], -2-2j, 2+2j, 1024, 1024, 1e-6, 64)

Sellest pildist saame aru, kui me vaatame $f'(z) = 5z^4 - 1$. Nullkohad on hulgas

$$C_0 = \left\{\frac{1}{\sqrt[4]{5}}, \frac{i}{\sqrt[4]{5}}, -\frac{1}{\sqrt[4]{5}}, -\frac{i}{\sqrt[4]{5}}\right\}.$$

Kui me vaatame näiteks $z_0 = \frac{i}{\sqrt[4]{5}}$ piirkonda, näeme et see punkt asub kõige nelja basseini puutumispunktis:

In [12]:
newtpoly([0, 1, -1, 1j, -1j], -0.128+0.54j, 0.128+0.7j, 1024, 640, 1e-6, 64)

Huvitav nähtus ilmub ka kui me valime näiteks $f(z) = z^3 -2z + 2$:

In [13]:
newtpoly([-1.7629, 0.8846+0.5897j, 0.8846-0.5897j], -3.0-2.5j, 2.0+2.5j, 1024, 1024, 1e-6, 64)

Siin on näha et näiteks $z = 0$ ja $z = 1$ ümbrused on musta värviga tähistatud, seega Newtoni meetod ei koondu, kui altingimus on valitud nendest hulkadest. Tõepoolest, kui me vaatame

$$F(z) = z - \frac{f(z)}{f'(z)} = -2\frac{z^3 - 1}{3z^2 - 2},$$

siis leiame, et

$$F(0) = 1, \quad F(1) = 0,$$

seega on olemas perioodiline trajektoor. Lisaks leiame, et

$$\left.\frac{d}{dz}F^2(z)\right|_{z = 0} = \left.\frac{d}{dz}F^2(z)\right|_{z = 1} = 0,$$

millest järeldub et ta on (isegi ülistabiilne) atraktor. Must värv näitab tema atraktiivsuse basseini. Siin on sellest ka suurem pilt:

In [14]:
newtpoly([-1.7629, 0.8846+0.5897j, 0.8846-0.5897j], -0.2-0.4j, 0.2+0.4j, 1024, 2048, 1e-6, 64)

Sarnane struktuur tekib, kui me valime näiteks

$$f(z) = \left(z - \frac{5}{2}\right)\left(z + \frac{5}{2}\right)(z - i)(z + i) = \left(z^2 - \frac{25}{4}\right)(z^2 + 1).$$

Pilt on selline:

In [15]:
newtpoly([2.5, -2.5, 1j, -1j], -3.2-2.0j, 3.2+2.0j, 1024, 640, 1e-6, 64)

Sellisel juhul on isegi mitu perioodilist trajektoori:

In [16]:
newtpoly([2.5, -2.5, 1j, -1j], -0.625j, 2+0.625j, 1024, 640, 1e-6, 64)

Trigonomeetrilised funktsioonid

Veel üks huvitav valik on $f(z) = \cos z$. Nullkohad on $A_0 = \pi\left(\mathbb{Z} + \frac{1}{2}\right)$. Basseini ääred on samuti fraktaalsed:

In [17]:
newton("cos(z)", [np.pi * (k + 0.5) for k in range (-10, 10)], -0.4-0.25j, 0.4+0.25j, 1024, 640, 1e-6, 256)

<lambdifygenerated-8>:2: RuntimeWarning: overflow encountered in sin
  return (-sin(z))
<lambdifygenerated-5>:2: RuntimeWarning: overflow encountered in sin
  return (-sin(z))