%%html
<style type="text/css">
input.good:checked + label {color: green}
input.bad:checked + label {color: red}
input.good:checked + label::after {color: green; content: ' Õige vastus!'}
input.bad:checked + label::after {color: red; content: ' Vale vastus!'}

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

Kaose klassikalised näited

Kõige tuntumad mittelineaarsed süsteemid, millel on kaootiline käitumine, on Lorenzi ja Rössleri süsteemid. Mõlemad on pidevad süsteemid kolmemõõtmelises faasiruumid.

Kasulikud Pythoni funktsioonid

Selleks et süsteemide trajektoore arvutada ja atraktoreid joonistada kasutame SciPy, NumPy ja PyPlot.

%matplotlib inline

from scipy.integrate import odeint
import matplotlib.pyplot as plt
import numpy as np
from mpl_toolkits.mplot3d import Axes3D

plt.rc('figure', figsize=(10.0, 7.5))

Lisaks defineerime funktsiooni, mis lahendab dünaamilise süsteemi võrrandeid ajavahemikus $[0, t_1]$ ja joonistab trajektoori kolmemõõtmelises ruumis ajavahemikus $[t_0, t_1]$:

def attrplot(func, par, ini, t0, t1, dt):
    n0 = (int)(t0 / dt)
    n1 = (int)(t1 / dt)
    T = np.linspace(0, t1, n1 + 1)
    sol = odeint(func, ini, T, args = par)
    xsol = sol[n0:, 0]
    ysol = sol[n0:, 1]
    zsol = sol[n0:, 2]
    fig = plt.figure()
    ax = fig.add_subplot(111, projection = '3d')
    ax.plot(xsol, ysol, zsol)
    plt.show()

Lorenzi süsteem

Definitsioon

Lorenzi süsteemi kirjeldavad võrrandid

\begin{align} \dot{x} &= \sigma(y - x),\\ \dot{y} &= x(\rho - z) - y,\\ \dot{z} &= xy - \beta z, \end{align}

kus $\sigma, \rho, \beta$ on konstantsed parameetrid. Seega defineerime:

def lorenz(X, t, s, b, r):
    x, y, z = X
    dx = s * (y - x)
    dy = r * x - y - x * z
    dz = x * y - b * z
    return [dx, dy, dz]

Tavaliselt valitakse $\sigma = 10$, $\beta = \frac{8}{3}$ ja $\rho$ varieeritakse. Siin on mõned kasulikud väärtused:

sigma = 10.0
beta = 8.0 / 3.0

rhos = [100.0, 99.8, 99.58, 99.0, 35.0, 28.0]

Simulatsioon

Algtingimuseks valime $x_0 = z_0 = 0$ ja $y_0 = 1$:

X0 = [0.0, 1.0, 0.0]

Järgmisena valime ajavahemiku. Simuleeritud ajavahemik on $[0, 5000]$. Sellest joonistame rajektoorina ainult osa $[4900, 5000]$, ja Poincaré lõike jaoks valime $[4500, 5000]$. Igat ajaühikut jagame 1000 arvutussammuks simulatsiooni jaoks.

tplot = 100
tvis = 500
tsim = 5000
steps = 1000

nsim = steps * tsim
nvis = steps * tvis
nplot = steps * tplot
t = np.linspace(0, tsim, nsim + 1)

Lõpuks simuleerime dünaamikat.

sols = [odeint(lorenz, X0, t, args = (sigma, beta, rho)) for rho in rhos]

Ruumiline joonis

Nüüd joonistame atraktori erinevate $\rho$ väärtuste puhul:

for sol in sols:
    fig = plt.figure()
    ax = fig.add_subplot(111, projection = '3d')
    ax.plot(sol[-nplot:, 0], sol[-nplot:, 1], sol[-nplot:, 2])
    plt.show()

Poincaré lõige ja kujutus

Järgmisena uurime Poincaré lõikeid. Selleks nihutame trajektoori mis vastab parameetri väärtusele $\rho = 28$ ühe ajasammu $\delta t$ võrra.

sol = sols[-1]
sola = sol[-nvis:-2, :]
solb = sol[-nvis + 1:-1, :]
solc = sol[-nvis + 2:, :]

Kõige lihtsam valik on uurida $z$ maksimumväärtusi, kus kehtib $xy = \beta z$. Numbriliselt leiame neid tingimuse kaudu, et $z(t)$ on suurem kui naabriväärtused $z(t - \delta t)$ ja $z(t + \delta t)$. Seejärel joonistame, kuidas järgmine $z$ maksimum sõltub eelmisest maksimumist. Tulemus meenutab telgikujutust:

cut = np.logical_and(sola[:, 2] < solb[:, 2], solb[:, 2] > solc[:, 2])
maxz = np.extract(cut, solb[:, 2])
zz = np.linspace(min(maxz), max(maxz), 101)
plt.plot(zz, zz, 'r')
plt.plot(maxz[:-1], maxz[1:], 'b.', markersize=5)
plt.grid()
plt.show()

Lisaks võime ka uurida kuju keskmist osa, kus $z = 31$ ja $\dot{z} < 0$. Sellisel juhul meenutab Poincaré kujutus saehambakujutust.

cutz = 31
cut = np.logical_and(np.logical_and(abs(solb[:, 2] - cutz) < abs(sola[:, 2] - cutz), abs(solb[:, 2] - cutz) < abs(solc[:, 2] - cutz)), sola[:, 2] > solc[:, 2])
extx = np.extract(cut, solb[:, 0])
xx = np.linspace(min(extx), max(extx), 101)

plt.plot(xx, xx, 'r')
plt.plot(extx[:-1], extx[1:], 'b.', markersize=5)
plt.grid()
plt.show()

Feigenbaumi diagramm

Feigenbaumi diagrammi joonistamiseks kasutame jälle $z$ koordinaadi maximumväärtusi.

def lorfeig(rhomin, rhomax, nrho):
    rhos = np.linspace(rhomin, rhomax, nrho + 1)
    vals = np.empty((0, 2))

    for rho in rhos:
        sol = odeint(lorenz, X0, t, args = (sigma, beta, rho))
        sola = sol[-nvis:-2, :]
        solb = sol[-nvis + 1:-1, :]
        solc = sol[-nvis + 2:, :]
        cut = np.logical_and(sola[:, 2] < solb[:, 2], solb[:, 2] > solc[:, 2])
        maxz = np.extract(cut, solb[:, 2])
        maxzz = np.transpose(np.array([np.full_like(maxz, rho), maxz]))
        vals = np.concatenate((vals, maxzz))

    plt.scatter(vals[:, 0], vals[:, 1], marker=".", s=1, lw=0)
    plt.xlabel("$\\rho$")
    plt.ylabel("z")
    plt.xlim(rhomin, rhomax)
    plt.grid()
    plt.show()

Selle funktsiooni abil joonistame Feigenbaumi diagrammi parameetri vahemikus $\rho \in [99.1, 100]$.

lorfeig(99.1, 100.0, 1500)

Edasi uurides leiame, et kaose vahel ilmuvad korrapärased vahemikud.

lorfeig(25.0, 100.0, 1500)

Rössleri süsteem

Definitsioon

Rössleri süsteemi kirjeldavad võrrandid

\begin{align} \dot{x} &= -y - z,\\ \dot{y} &= x + ay,\\ \dot{z} &= b + z(x - c), \end{align}

kus $a, b, c$ on konstantsed parameetrid. Seega defineerime:

def roessler(X, t, a, b, c):
    x, y, z = X
    dx = -y - z
    dy = x + a * y
    dz = b + x * z - c * z
    return [dx, dy, dz]

Tavaliselt valitakse $a = b = \frac{1}{5}$ ja $c$ varieeritakse.

aa = 0.2
bb = 0.2

ccs = [3.5, 4.0, 4.2, 5.0, 5.3, 5.7, 6.0, 6.7, 7.5, 8.0]

Simulatsioon

Algtingimuseks valime $y_0 = z_0 = 0$ ja $x_0 = -1$:

X0 = [-1.0, 0.0, 0.0]

Järgmisena valime ajavahemiku. Simuleeritud ajavahemik on $[0, 5000]$. Sellest joonistame rajektoorina ainult osa $[4900, 5000]$, ja Poincaré lõike jaoks valime $[4500, 5000]$. Igat ajaühikut jagame 1000 arvutussammuks simulatsiooni jaoks.

tplot = 500
tvis = 1000
tsim = 5000
steps = 1000

nsim = steps * tsim
nvis = steps * tvis
nplot = steps * tplot
t = np.linspace(0, tsim, nsim + 1)

Lõpuks simuleerime dünaamikat.

sols = [odeint(roessler, X0, t, args = (aa, bb, cc)) for cc in ccs]

Ruumiline joonis

Nüüd joonistame atraktori erinevate $c$ väärtuste puhul:

for sol in sols:
    fig = plt.figure()
    ax = fig.add_subplot(111, projection = '3d')
    ax.plot(sol[-nplot:, 0], sol[-nplot:, 1], sol[-nplot:, 2])
    plt.show()

Poincaré kujutus

Ühe $c$ väärtuse jaoks näitame ka Poincaré lõiget. Selleks valime $x$ koordinaadi miinimumväärtusi.

sol = sols[5]
sola = sol[-nvis:-2, :]
solb = sol[-nvis + 1:-1, :]
solc = sol[-nvis + 2:, :]

Saadud Poincarè kujutus meenutab logistilist kujutust.

cut = np.logical_and(sola[:, 0] > solb[:, 0], solb[:, 0] < solc[:, 0])
minx = -np.extract(cut, solb[:, 0])
xx = np.linspace(min(minx), max(minx), 101)
plt.plot(xx, xx, 'r')
plt.plot(minx[:-1], minx[1:], 'b.', markersize=5)
plt.grid()
plt.show()

Feigenbaumi diagramm

Lõpuks joonistame ka Feigenbaumi diagrammi parameetri vahemikus $c \in [2.5, 10.0]$.

ccmin = 2.5
ccmax = 10.0
ncc = 1500

ccs = np.linspace(ccmin, ccmax, ncc + 1)

Feigenbaumi diagrammi joonistamiseks kasutame jälle $x$ koordinaadi miinimumväärtusi.

vals = np.empty((0, 2))

for cc in ccs:
    sol = odeint(roessler, X0, t, args = (aa, bb, cc))
    sola = sol[-nvis:-2, :]
    solb = sol[-nvis + 1:-1, :]
    solc = sol[-nvis + 2:, :]
    cut = np.logical_and(sola[:, 0] > solb[:, 0], solb[:, 0] < solc[:, 0])
    minx = -np.extract(cut, solb[:, 0])
    minxx = np.transpose(np.array([np.full_like(minx, cc), minx]))
    vals = np.concatenate((vals, minxx))

plt.scatter(vals[:, 0], vals[:, 1], marker=".", s=1, lw=0)
plt.xlabel("c")
plt.ylabel("-x")
plt.xlim(ccmin, ccmax)
plt.grid()
plt.show()