In [1]:
%%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!'}
</style>

Poincaré kujutus ja lõige

Mõisted

Kui $n$-mõõtmelise dünaamilise süsteemi faasiruumil $Q$ on selline $(n - 1)$-mõõtmeline alamhulk $S \subset Q$, et igale algtingimusele $\underline{x}_0 \in S$ vastav trajektoor $t \mapsto \underline{x}(t)$, kus $\underline{x}(0) = \underline{x}_0$, läbib hulka $S$ uuesti ajahetkel $T > 0$, niiet $\underline{x}(T) \in S$, ja dünaamikast $f$ määratud tuletis

$$\dot{\underline{x}}(0) = \left.\frac{d}{dt}\underline{x}(t)\right|_{t = 0} = f(\underline{x}_0)$$

aja järgi on $S$ suhtes transversaalne (mitte paralleelne), siis nimetatakse alamhulka $S$ Poincaré lõikeks.

Tingimusest, et $\dot{\underline{x}}(0)$ ja $S$ on transversaalsed, järeldub et punkti $\underline{x}_0 \in S$ läbiv trajektoor ei ole $S$ suhtes paralleelne. Seega võib valida aega $T$ nii, et $\underline{x}(0) \in S$, $\underline{x}(T) \in S$ ja $\underline{x}(t) \notin S$ kui $0 < t < T$. Kui $T$ on sellisena valitud, siis nimetatakse funktsiooni

$$\begin{array}{rcrcl} P & : & S & \to & S\\ && \underline{x}_0 & \mapsto & \underline{x}(T) \end{array}$$

esimese tagasituleku kujutuseks või Poincaré kujutuseks.

Üleminek diskreetsetele süsteemidele

Kujutus $P: S \to S$ määrab diskreetse dünaamilise süsteemi dünaamikat. Süsteemi faasiruum on $S$, ja trajektoor, mis vastab algtingimusele $\underline{x}_0 \in S$, on jada $(\underline{x}_n, n \in \mathbb{N})$, mis rahuldab $\underline{x}_{n + 1} = P(\underline{x}_n)$. Pidev süsteem ja Poincaré lõike abil tuletatud diskreetne süsteem on tihedalt omavahel seotud:

  • Kui diskreetsel süsteemil on püsipunkt $\underline{x}_* \in S$, mis rahuldab $P(\underline{x}_*) = \underline{x}_*$, siis on pideva süsteemi trajektoor, mis vastab algtingimusele $\underline{x}_0 = \underline{x}_*$, perioodiline. Trajektoori periood on aeg $T$ mis mõõdub kuni trajektoor taas lõikab hulka $S$.
  • Kui diskreetse süsteemi trajektoor on perioodiline, see tähendab et on olemas $N \in \mathbb{N}$ nii, et kehtib $\underline{x}_N = \underline{x}_0$, siis on pideva süsteemi trajektoor, mis vastab algtingimusele $\underline{x}_0$, samuti perioodiline.
  • Kui diskreetse süsteemi trajektoor ei ole perioodiline, siis pideva süsteemi trajektoor ka ei ole perioodiline.

Selleks et uurida pideva süsteemi trajektoore, nende perioodilisust ja stabiilsust, on tihti kasulik uurida hoopis diskreetset süsteemi. Sellepärast tegeleme järgnevates loengutes diskreetsete süsteemidega.

Interaktiivne näide: hõõrdejõuga mittelineaarne pendel ja erinevad Poincaré lõiked

Järgnevalt on selgitatud, kuidas on võimalik leida Poincaré lõiget $S$ ja konstrueerida Poincaré kujutust $P$. Selleks kasutame SciPy, NumPy ja PyPlot.

In [2]:
%matplotlib inline

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

Süsteemi definitsioon

Süsteem, mida me uurime, on hõõrdejõuga füüsikaline pendel, mille tagasitõmmejõud on mittelineaarne. Pendli asukoht $x \in S^1$ on punkt ringjoonel, kiirus $v \in \mathbb{R}$ on reaalarv. Faasiruum on seega $Q = S^1 \times \mathbb{R}$. Pendli liikumist kirjeldab võrrandite süsteem:

  • $\dot{x} = v,$
  • $\dot{v} = -\frac{hv + k\sin x}{m},$

kus $m$ on pendli mass, $k$ on vedru parameeter ja $h$ on hõõrdetegur. Järgnev Pythoni funktsioon kirjeldab süsteemi dünaamikat:

In [3]:
def pendel(y, t, h, k, m):
    x, v = y
    dydt = [v, -(h * v + k * np.sin(x)) / m]
    return dydt

Konstandid, mida me kasutame, on näiteks sellised:

In [4]:
h = 0.5
k = 0.5
m = 1

Faasiportree

Järgmisena joonistame faasiportreet. Seda saab teha PyPlot funktsioonidega streamplot ja quiver. Paneme tähele, et $x \in S^1$, ja vastav koordinaat asub näiteks vahemikus $(-\pi, \pi]$. Tulemus on selline:

In [5]:
fig = plt.figure()
ax = fig.add_subplot(111)
ax.set_aspect('equal', 'box')
ax.set_xlabel('x')
ax.set_ylabel('v')

X = np.linspace(-np.pi, np.pi, 101)
Y = np.linspace(-3.5, 3.5, 101)
XX, YY = np.meshgrid(X, Y)
DX, DY = pendel([XX, YY], 0, h, k, m)
ax.streamplot(XX, YY, DX, DY, color='b')

X = np.linspace(-np.pi, np.pi, 21)
Y = np.linspace(-3.5, 3.5, 21)
XX, YY = np.meshgrid(X, Y)
DX, DY = pendel([XX, YY], 0, h, k, m)
ax.quiver(XX, YY, DX, DY, color='r')

plt.show()

Faasiportree näitab, et on kolm trajektooritüüpi:

  • Tasakaalupunktid $(0, 0)$ ja $(\pi, 0)$ on püsipunktid. Nendele algtingimustele vastavad trajektoorid koosnevad ainult ühest punktist. Alumine püsipunkt $(0, 0)$ on atraktor, ülemine püsipunkt $(\pi, 0)$ on sadul.
  • Teised trajektoorid on üldjuhul sarnased: Pendel ületab tippu $x = \pm\pi$ nii kaua kui energiast piisab, ja siis võngub püsipunkti $(0, 0)$ lähenedes.
  • Erijuht on kaks trajektoori, mis lähenevad sadulale $(\pi, 0)$ kas vasakult või paremalt.

Ehijuhi trajektoore on mõistlik ka joonistada, pildil on nad punased:

In [6]:
# Kõige lihtsam on neid trajektoore leida, kui me sadulast tagasi integreerime.
def rpendel(y, t, h, k, m):
    return [-dy for dy in pendel(y, t, h, k, m)]

fig = plt.figure()
ax = fig.add_subplot(111)
ax.set_aspect(0.5)
ax.set_xlabel('x')
ax.set_ylabel('v')

# Joonistame kõigepealt süsteemi voolu.
X = np.linspace(-np.pi, np.pi, 101)
Y = np.linspace(-7, 7, 201)
XX, YY = np.meshgrid(X, Y)
DX, DY = pendel([XX, YY], 0, h, k, m)
ax.streamplot(XX, YY, DX, DY, color='b')

# Nüüd integreerime liikumisvõrrandeid.
t = np.linspace(0, 16, 1000)
ysolm = odeint(rpendel, [np.pi - 0.00001, 0], t, args = (h, k, m))
ysolp = odeint(rpendel, [np.pi + 0.00001, 0], t, args = (h, k, m))

# Pendli asukoha lükkame vahemikku [-pi, pi].
xsolm = np.mod(ysolm[:, 0] + np.pi, 2 * np.pi) - np.pi
xsolp = np.mod(ysolp[:, 0] + np.pi, 2 * np.pi) - np.pi

# Jooonistamiseks lõikame trajektoori ka tükkideks, sest meil on lõigatud faasiruum.
cutm = np.argwhere(np.fabs(xsolm[1:] - xsolm[:-1]) > np.pi) + 1
cutp = np.argwhere(np.fabs(xsolp[1:] - xsolp[:-1]) > np.pi) + 1

ycutm = np.split(ysolm, cutm[:, 0])
ycutp = np.split(ysolp, cutp[:, 0])
ycut = ycutm + ycutp

# Lõpuks joonistame.
for yy in ycut:
    xx = np.mod(yy[:, 0] + np.pi, 2 * np.pi) - np.pi
    vv = yy[:, 1]
    ax.plot(xx, vv, 'r')

plt.show()

Poincaré lõiked

Poincaré lõike mõiste nõuab, et trajektooritel, mis vastavad algtingimusele $\underline{x}(0) = \underline{x}_0 \in S$, on transversaalne puutevektor $\dot{\underline{x}}(0)$. Seega püsipunktid ei saa kuuluda hulka $S$.

Eelnevast järeldub, et meil on mitu võimalust valida Poincaré lõiget. Seda saame teha, kui me vaatame süsteemi trajektoore. Selleks on kasulik defineerida sellist funktsiooni, mis joonistab trajektoori ja Poincaré lõiget, ja selle abil uurida erinevaid trajektoore ja algtingimusi. Trajektoor, mida me uurime, on roheline, ja valitud Poincaré lõige on pruun.

In [7]:
def tplot(y0, pl):
    fig = plt.figure()
    ax = fig.add_subplot(111)
    ax.set_aspect(0.5)
    ax.set_xlabel('x')
    ax.set_ylabel('v')

    ax.streamplot(XX, YY, DX, DY, color='b')

    for po in np.array(pl):
        ax.plot(po[:, 0], po[:, 1], '#804000')

    for yy in ycut:
        xx = np.mod(yy[:, 0] + np.pi, 2 * np.pi) - np.pi
        vv = yy[:, 1]
        ax.plot(xx, vv, 'r')

    t = np.linspace(0, 20, 1000)
    sol = odeint(pendel, y0, t, args = (h, k, m))
    xsol = np.mod(sol[:, 0] + np.pi, 2 * np.pi) - np.pi
    cut = np.argwhere(np.fabs(xsol[1:] - xsol[:-1]) > np.pi) + 1
    cuts = np.split(sol, cut[:, 0])

    for yy in cuts:
        xx = np.mod(yy[:, 0] + np.pi, 2 * np.pi) - np.pi
        vv = yy[:, 1]
        ax.plot(xx, vv, '#008000')

    plt.show()

1. $x = 0, 0 < |v| < v_1$

Kui me valime algtingimust $(0, v)$ nii, et $0 < |v| < v_1$, kus $v_1$ on see kiirus, et pendel jõuab täpselt sadulani $(\pi, 0)$, siis pendel võngub tagasi, ületab alumist punkti $x = 0$ teises suunas väiksema kiirusega, ja siis veel kord vahetab suunda ja ületab punkti $x = 0$ jälle algses suunas veel väiksema kiirusega $v'$:

In [8]:
tplot([0, 2], [])

1a.

$$S^{v+}_1 = \{(0, v), 0 < v < v_1\}$$

Trajektoor, mis alustab sellest hulgast $S^{v+}_1$, jõuab hulka $S^{v+}_1$ tagasi, kui pendel on kaks korda oma suunda vahetanud ja liigub jälle algses suunas läbi punkti $x = 0$:

In [9]:
tplot([0, 2], [[(0, 0), (0, 2.5)]])

1b.

$$S^{v-}_1 = \{(0, v), -v_1 < v < 0\}$$

Siin kehtib täpselt sama kui eelmise valiku puhul, lihtsalt et algne kiirus on negatiivne:

In [10]:
tplot([0, -2], [[(0, 0), (0, -2.5)]])

1c.

$$S^{v\pm}_1 = \{(0, v), 0 < |v| < v_1\}$$

Trajektoor, mis algab sellest hulgast $S^{v\pm}_1$, jõuab hulka $S^{v\pm}_1$ tagasi, kui pendel on ükskord oma suunda vahetanud ja liigub vastassuunas läbi punkti $x = 0$.

In [11]:
tplot([0, 2], [[(0, -2.5), (0, 2.5)]])

2. $x = 0, v_n < |v| < v_{n + 1}$

Lisaks kiirustele $v_1$ ja $-v_1$ viivad sadulani ka algtingimused $(0, v_n)$, kus $n \in \mathbb{N}$, ja pendel ületab tippu $x = \pi$ veel $n - 1$ korda enne kui ta läheneb sadulale. Teistele algtingimustele kujul $(0, v)$ vastavad sellised trajektoorid:

  • Kui $0 < |v| < v_1$, siis pendel lihtsalt võngub - seda me juba uurisime.
  • Kui $v_n < |v| < v_{n + 1}$, siis pendel ületab tippu $n$ korda samas suunas, enne kui ta jääb võnkuma.

Teine trajektoori tüüp on siis selline:

In [12]:
tplot([0, 5], [])

2a.

$$S^{v+}_n = \{(0, v), 0 < v < v_n\} \setminus \{(0, v_1), \ldots, (0, v_{n - 1})\}$$

Trajektoor, mis algab sellest hulgast $S^{v+}_n$ kiirusega $v$, mis rahuldab $v_k < v < v_{k + 1}$, kus $k < n$, jõuab hulka $S^{v+}_n$ tagasi, kui pendel on ükskord tippu ületanud ja liigub väiksema kiirusega edasi algses suunas läbi punkti $x = 0$. Kui algne kiirus $v$ rahuldab hoopis $0 < v < v_1$, kehtib sama kui seletatud punkti 1a all. Alljärgnev joonis näitab $S^{v+}_2$.

In [13]:
tplot([0, 5], [[(0, 0), (0, 5.55)]])

2a'.

$$S^{v+}_{\infty} = \bigcup_{n = 1}^{\infty}S^{v+}_n = \{(0, v), v > 0\} \setminus \{(0, v_n), n \in \mathbb{N}\}$$

Poincaré lõikeks sobib ka terve positiivne poolsirge, millest on eemaldatud need punktid, kust trajektoor läheneb sadulale.

In [14]:
tplot([0, 7], [[(0, 0), (0, 8)]])

2b.

$$S^{v-}_n = \{(0, v), -v_n < v < 0\} \setminus \{(0, -v_1), \ldots, (0, -v_{n - 1})\}$$

Siin kehtib kvalitatiivselt sama kui hulga $S^{v+}_n$ puhul, lihtsalt et algne kiirus on negatiivne. Alljärgnev joonis näitab $S^{v-}_2$.

In [15]:
tplot([0, -5], [[(0, 0), (0, -5.55)]])

2b'.

$$S^{v-}_{\infty} = \bigcup_{n = 1}^{\infty}S^{v-}_n = \{(0, v), v < 0\} \setminus \{(0, -v_n), n \in \mathbb{N}\}$$

Poincaré lõikeks sobib ka terve negatiivne poolsirge, millest on eemaldatud need punktid, kust trajektoor läheneb sadulale.

In [16]:
tplot([0, -7], [[(0, 0), (0, -8)]])

2c.

$$S^{v\pm}_n = \{(0, v), 0 < |v| < v_n\} \setminus \{(0, \pm v_1), \ldots, (0, \pm v_{n - 1})\}$$

Trajektoor, mis algab sellest hulgast $S^{v\pm}_n$ kiirusega $v$, mis rahuldab $v_k < v < v_{k + 1}$, kus $k < n$, jõuab hulka $S^{v\pm}_n$ tagasi, kui pendel on ükskord tippu ületanud ja liigub väiksema kiirusega edasi algses suunas läbi punkti $x = 0$. Kui algne kiirus $v$ rahuldab hoopis $0 < |v| < v_1$, kehtib sama kui seletatud punkti 1c all. Alljärgnev joonis näitab $S^{v\pm}_2$.

In [17]:
tplot([0, 5], [[(0, -5.55), (0, 5.55)]])

2c'.

$$S^{v\pm}_{\infty} = \bigcup_{n = 1}^{\infty}S^{v\pm}_n = \{(0, v), v \neq 0\} \setminus \{(0, \pm v_n), n \in \mathbb{N}\}$$

Poincaré lõikeks sobib ka terve sirge $x = 0$, millest on eemaldatud need punktid, kust trajektoor läheneb sadulale.

In [18]:
tplot([0, 7], [[(0, -8), (0, 8)]])

3. $0 < |x| < \pi, v = 0$

Kui me valime algtingimust $(x, 0)$ nii, et pendel on paigalolekus, aga ei asu püsipunktis, siis ta hakkab võnkuma, ületab alumist punkti $x = 0$, vahetab suunda, ületab punkti $x = 0$ vastassuunas jne:

In [19]:
tplot([3, 0], [])

3a.

$$S^{x+} = \{(x, 0), 0 < x < \pi\}$$

Trajektoor, mis algab sellest hulgast $S^{x+}$, jõuab sinna tagasi, kui pendel on võnkunud ühes suunas $v < 0$ ja tagasi teises suunas $v > 0$, kuni ta on jälle paigalolekus ja vahetab teist korda liikumissuunda.

In [20]:
tplot([3, 0], [[(0, 0), (np.pi, 0)]])

3b.

$$S^{x-} = \{(x, 0), -\pi < x < 0\}$$

Trajektoor, mis algab sellest hulgast $S^{x-}$, jõuab sinna tagasi, kui pendel on võnkunud ühes suunas $v > 0$ ja tagasi teises suunas $v < 0$, kuni ta on jälle paigalolekus ja vahetab teist korda liikumissuunda.

In [21]:
tplot([-3, 0], [[(0, 0), (-np.pi, 0)]])

3c.

$$S^{x\pm} = \{(x, 0), 0 < |x| < \pi\}$$

Veel üks võimalus on valida hulka $S^{x\pm}$ nii, et pendel on paigalolekus, ja asub kas vasakul või paremal pool. Trajektoor, mis algab sellest hulgast, jõuab sinna tagasi pärast esimest võnkumist, kui ta on jälle paigalolekus, kuid siis teisel pool tasakaalupunktist.

In [22]:
tplot([3, 0], [[(-np.pi, 0), (np.pi, 0)]])

Poincaré kujutus

Lõpuks konstrueerime Poincaré kujutusi erinevate Poincaré lõigete valiku puhul. Täpsemalt uurime valikuid $S^{v+}_{\infty}$, $S^{v\pm}_{\infty}$, $S^{x+}$ ja $S^{x\pm}$.

$S^{v+}_{\infty}$

Kõigepealt uurime Poincaré lõiget

$$S^{v+}_{\infty} = \{(0, v), v > 0\} \setminus \{(0, v_n), n \in \mathbb{N}\},$$

ja konstrueerime Poincaré kujutust. Poincaré kujutuse $P$ kvalitativset käitumist saame leida, kui me vaatame erinevaid trajektoore. Sõltuvalt algtingimusest $(0, v)$ leiame:

  • Kui $0 < v < v_1$, siis pendel jõuab tagasi hulka $S^{v+}_{\infty}$ siis, kui ta on edasi ja tagasi võnkunud, ja liigub taas algses suunas, kuid väiksemal kiirusel $P(v)$, mis rahuldab $0 < P(v) < v$.
  • Kui $v_1 < v < v_2$, siis pendel ületab tippu ja jõuab tagasi hulka $S^{v+}_{\infty}$ kiirusega $P(v)$, mis on tipu ületamiseks liiga väike, ehk $0 < P(v) < v_1$.
  • Kui $v_n < v < v_{n + 1}$ ja $n > 1$, siis pendel ületab tippu ja jõuab tagasi hulka $S^{v+}_{\infty}$ kiirusega $P(v)$, millega ta saab tippu veel $n - 1$ korda ületada, ehk $v_{n - 1} < P(v) < v_n$.

Täpsemalt saame uurida Poincaré kujutust järgneva funktsiooni abil.

In [23]:
def poinca(v0):
    y0 = [0, v0]

    # Eeldame, et trajektoor jõuab tagasi samasse hulka 20s jooksul, ja integreerime.
    t = np.linspace(0, 20, 1000)
    sol = odeint(pendel, y0, t, args = (h, k, m))
    
    # Jagame lahendust x ja v väärtusteks.
    xsol = np.mod(sol[:, 0] + np.pi, 2 * np.pi) - np.pi
    vsol = sol[:, 1]
    
    # Otsime, millal:
    # 1. Asukoha x absoluutväärtus on minimaalne (pendel ületab alumist punkti).
    # 2. Kiiruse v märk on sama kui algsel hetkel.
    xabs = np.fabs(xsol)
    xmin = np.r_[True, xabs[1:] < xabs[:-1]] & np.r_[xabs[:-1] < xabs[1:], True] & (vsol * v0 > 0)
    xzero = np.argwhere(xmin)
    first = xzero[1, 0]
    t0 = t[first]

    # Arvutame veel kord sobitatud ajasammuga, et tulemus oleks täpsem.
    t = np.linspace(0, 1.5 * t0, 1000)
    sol = odeint(pendel, y0, t, args = (h, k, m))
    xsol = np.mod(sol[:, 0] + np.pi, 2 * np.pi) - np.pi
    vsol = sol[:, 1]
    xabs = np.fabs(xsol)
    xmin = np.r_[True, xabs[1:] < xabs[:-1]] & np.r_[xabs[:-1] < xabs[1:], True] & (vsol * v0 > 0)
    xzero = np.argwhere(xmin)
    first = xzero[1, 0]

    # Tulemus on tagasituleku aeg T(v) ja Poincaré kujutus P(v).
    return [t[first], vsol[first]]

Selle funktsiooni abil võime nüüd $T(v)$ ja $P(v)$ arvutada.

In [24]:
vv = np.linspace(0.01, 10, 500)
res = np.array([poinca(v) for v in vv])
tt = res[:, 0]
pp = res[:, 1]

Joonistame kõigepealt $T(v)$. Joonis näitab, et kiiruse $v = v_1$ ümbruses tagasitulekuaeg hajub, $T(v) \to \infty$. See on sellepärast, et sadula $(\pi, 0)$ ümbruses liigub pendel aeglaselt, ja võnkumine võtab kauem aega.

In [25]:
fig = plt.figure()
ax = fig.add_subplot(111)
ax.set_xlabel('v')
ax.set_ylabel('T(v)')
ax.plot(vv, tt)
ax.grid()
plt.show()

Järgmisena joonistame $P(v)$. Paneme tähele, et kehtib $0 < P(v) < v$, kuna pendel aeglustub, ja Poincaré lõige on valitud sellisena, et liikussuund on alati positiivne. Joonis näitab, et $P(v)$ ei ole pidev kiiruse $v = v_1$ ümbruses. Aga tuletame meelde, et vastav algtingimus $(0, v_1) \notin S^{v+}_{\infty}$. Seega on Poincaré kujutus siiski pidev oma lähtehulgas $S^{v+}_{\infty}$.

In [26]:
fig = plt.figure()
ax = fig.add_subplot(111)
ax.set_xlabel('v')
ax.set_ylabel('P(v)')
ax.plot(vv, pp)
ax.grid()
plt.show()

$S^{v\pm}_{\infty}$

Järgmisena uurime Poincaré lõiget

$$S^{v\pm}_{\infty} = \{(0, v), v \neq 0\} \setminus \{(0, \pm v_n), n \in \mathbb{N}\},$$

ja konstrueerime Poincaré kujutust. Poincaré kujutuse $P$ kvalitativset käitumist saame leida, kui me vaatame erinevaid trajektoore. Sõltuvalt algtingimusest $(0, v)$ leiame:

  • Kui $0 < |v| < v_1$, siis pendel jõuab tagasi hulka $S^{v\pm}_{\infty}$ siis, kui ta on oma liikumissuunda vahetanud, ja ületab punkti $x = 0$ kiirusel $P(v)$, mis rahuldab $0 < |P(v)| < v$ ja $\mathrm{sgn} P(v) = -\mathrm{sgn} v$.
  • Kui $v_1 < |v| < v_2$, siis pendel ületab tippu ja jõuab tagasi hulka $S^{v+}_{\infty}$ kiirusega $P(v)$, mis on tipu ületamiseks liiga väike, ehk $0 < |P(v)| < v_1$, aga suund jääb samaks, $\mathrm{sgn} P(v) = -\mathrm{sgn} v$.
  • Kui $v_n < |v| < v_{n + 1}$ ja $n > 1$, siis pendel ületab tippu ja jõuab tagasi hulka $S^{v+}_{\infty}$ kiirusega $P(v)$, millega ta saab tippu veel $n - 1$ korda ületada, ehk $v_{n - 1} < |P(v)| < v_n$, ja suund jääb samaks, $\mathrm{sgn} P(v) = -\mathrm{sgn} v$.

Täpsemalt saame uurida Poincaré kujutust järgneva funktsiooni abil.

In [27]:
def poincb(v0):
    y0 = [0, v0]

    # Eeldame, et trajektoor jõuab tagasi samasse hulka 20s jooksul, ja integreerime.
    t = np.linspace(0, 20, 1000)
    sol = odeint(pendel, y0, t, args = (h, k, m))
    
    # Jagame lahendust x ja v väärtusteks.
    xsol = np.mod(sol[:, 0] + np.pi, 2 * np.pi) - np.pi
    vsol = sol[:, 1]
    
    # Otsime, millal asukoha x absoluutväärtus on minimaalne (pendel ületab alumist punkti).
    xabs = np.fabs(xsol)
    xmin = np.r_[True, xabs[1:] < xabs[:-1]] & np.r_[xabs[:-1] < xabs[1:], True]
    xzero = np.argwhere(xmin)
    first = xzero[1, 0]
    t0 = t[first]

    # Arvutame veel kord sobitatud ajasammuga, et tulemus oleks täpsem.
    t = np.linspace(0, 1.5 * t0, 1000)
    sol = odeint(pendel, y0, t, args = (h, k, m))
    xsol = np.mod(sol[:, 0] + np.pi, 2 * np.pi) - np.pi
    vsol = sol[:, 1]
    xabs = np.fabs(xsol)
    xmin = np.r_[True, xabs[1:] < xabs[:-1]] & np.r_[xabs[:-1] < xabs[1:], True]
    xzero = np.argwhere(xmin)
    first = xzero[1, 0]

    # Tulemus on tagasituleku aeg T(v) ja Poincaré kujutus P(v).
    return [t[first], vsol[first]]

Selle funktsiooni abil võime nüüd $T(v)$ ja $P(v)$ arvutada.

In [28]:
vv = np.linspace(0.01, 10, 500)
vv = np.concatenate((np.flipud(-vv), vv), 0)
res = np.array([poincb(v) for v in vv])
tt = res[:, 0]
pp = res[:, 1]

Joonistame kõigepealt $T(v)$. Joonis näitab, et kiiruse $v = \pm v_1$ ümbruses tagasitulekuaeg hajub, $T(v) \to \infty$. See on sellepärast, et sadula $(\pi, 0)$ ümbruses liigub pendel aeglaselt, ja võnkumine võtab kauem aega.

In [29]:
fig = plt.figure()
ax = fig.add_subplot(111)
ax.set_xlabel('v')
ax.set_ylabel('T(v)')
ax.plot(vv, tt)
ax.grid()
plt.show()

Järgmisena joonistame $P(v)$. Paneme tähele, et kehtib $0 < |P(v)| < |v|$, kuna pendel aeglustub, ja vahetab suunda täpselt siis kui $0 < v < v_1$. Joonis näitab, et $P(v)$ ei ole pidev kiiruse $v = \pm v_1$ ümbruses. Aga tuletame veel kord meelde, et vastavad algtingimused $(0, \pm v_1) \notin S^{v\pm}_{\infty}$. Seega on Poincaré kujutus siiski pidev oma lähtehulgas $S^{v\pm}_{\infty}$.

In [30]:
fig = plt.figure()
ax = fig.add_subplot(111)
ax.set_xlabel('v')
ax.set_ylabel('P(v)')
ax.plot(vv, pp)
ax.grid()
plt.show()

$S^{x+}$

Uurime nüüd Poincaré lõiget

$$S^{x+} = \{(x, 0), 0 < x < \pi\},$$

ja konstrueerime Poincaré kujutust. Poincaré kujutuse $P$ kvalitativne käitumine on tegelikult lihtne:

  • Kui $0 < x < \pi$, siis pendel jõuab tagasi hulka $S^{x+}$ siis, kui ta on edasi ja tagasi võnkunud, ja vahetab suunda teist korda. Tema asukoht on siis $P(x)$, ja rahuldab $0 < P(x) < x$.

Täpsemalt saame uurida Poincaré kujutust järgneva funktsiooni abil.

In [31]:
def poincc(x0):
    y0 = [x0, 0]

    # Eeldame, et trajektoor jõuab tagasi samasse hulka 20s jooksul, ja integreerime.
    t = np.linspace(0, 20, 1000)
    sol = odeint(pendel, y0, t, args = (h, k, m))
    
    # Jagame lahendust x ja v väärtusteks.
    xsol = np.mod(sol[:, 0] + np.pi, 2 * np.pi) - np.pi
    vsol = sol[:, 1]
    
    # Otsime, millal:
    # 1. Kiiruse v absoluutväärtus on minimaalne (pendel vahetab suunda).
    # 2. Asukoha x märk on sama kui algsel hetkel.
    vabs = np.fabs(vsol)
    vmin = np.r_[True, vabs[1:] < vabs[:-1]] & np.r_[vabs[:-1] < vabs[1:], True] & (xsol * x0 > 0)
    vzero = np.argwhere(vmin)
    first = vzero[1, 0]
    t0 = t[first]

    # Arvutame veel kord sobitatud ajasammuga, et tulemus oleks täpsem.
    t = np.linspace(0, 1.5 * t0, 1000)
    sol = odeint(pendel, y0, t, args = (h, k, m))
    xsol = np.mod(sol[:, 0] + np.pi, 2 * np.pi) - np.pi
    vsol = sol[:, 1]
    vabs = np.fabs(vsol)
    vmin = np.r_[True, vabs[1:] < vabs[:-1]] & np.r_[vabs[:-1] < vabs[1:], True] & (xsol * x0 > 0)
    vzero = np.argwhere(vmin)
    first = vzero[1, 0]

    # Tulemus on tagasituleku aeg T(x) ja Poincaré kujutus P(x).
    return [t[first], xsol[first]]

Selle funktsiooni abil võime nüüd $T(x)$ ja $P(x)$ arvutada.

In [32]:
xx = np.linspace(0.01, np.pi - 0.01, 500)
res = np.array([poincc(x) for x in xx])
tt = res[:, 0]
pp = res[:, 1]

Joonistame kõigepealt $T(x)$. Joonis näitab, et tipu $x = \pi$ ümbruses tagasitulekuaeg hajub, $T(x) \to \infty$. See on sellepärast, et sadula $(\pi, 0)$ ümbruses liigub pendel aeglaselt, ja võnkumine võtab kauem aega.

In [33]:
fig = plt.figure()
ax = fig.add_subplot(111)
ax.set_xlabel('x')
ax.set_ylabel('T(x)')
ax.plot(xx, tt)
ax.grid()
plt.show()

Järgmisena joonistame $P(x)$. Paneme tähele, et kehtib $0 < P(x) < x$, kuna pendel kaotab energiat, ja Poincaré lõige on valitud sellisena, et pendli asukoht $x$ on alati positiivne.

In [34]:
fig = plt.figure()
ax = fig.add_subplot(111)
ax.set_xlabel('x')
ax.set_ylabel('P(x)')
ax.plot(xx, pp)
ax.grid()
plt.show()

$S^{x\pm}$

Lõpuks uurime Poincaré lõiget

$$S^{x\pm} = \{(x, 0), 0 < |x| < \pi\},$$

ja konstrueerime Poincaré kujutust. Poincaré kujutuse $P$ kvalitativne käitumine on ka sel juhul lihtne:

  • Kui $0 < |x| < \pi$, siis pendel jõuab tagasi hulka $S^{x\pm}$ siis, kui ta on teisele poole võnkunud, ja vahetab suunda esimest korda. Tema asukoht on siis $P(x)$, ja rahuldab $0 < |P(x)| < |x|$ ja $\mathrm{sgn} P(x) = -\mathrm{sgn} x$.

Täpsemalt saame uurida Poincaré kujutust järgneva funktsiooni abil.

In [35]:
def poincd(x0):
    y0 = [x0, 0]

    # Eeldame, et trajektoor jõuab tagasi samasse hulka 20s jooksul, ja integreerime.
    t = np.linspace(0, 20, 1000)
    sol = odeint(pendel, y0, t, args = (h, k, m))
    
    # Jagame lahendust x ja v väärtusteks.
    xsol = np.mod(sol[:, 0] + np.pi, 2 * np.pi) - np.pi
    vsol = sol[:, 1]
    
    # Otsime, millal kiiruse v absoluutväärtus on minimaalne (pendel vahetab suunda).
    vabs = np.fabs(vsol)
    vmin = np.r_[True, vabs[1:] < vabs[:-1]] & np.r_[vabs[:-1] < vabs[1:], True]
    vzero = np.argwhere(vmin)
    first = vzero[1, 0]
    t0 = t[first]

    # Arvutame veel kord sobitatud ajasammuga, et tulemus oleks täpsem.
    t = np.linspace(0, 1.5 * t0, 1000)
    sol = odeint(pendel, y0, t, args = (h, k, m))
    xsol = np.mod(sol[:, 0] + np.pi, 2 * np.pi) - np.pi
    vsol = sol[:, 1]
    vabs = np.fabs(vsol)
    vmin = np.r_[True, vabs[1:] < vabs[:-1]] & np.r_[vabs[:-1] < vabs[1:], True]
    vzero = np.argwhere(vmin)
    first = vzero[1, 0]

    # Tulemus on tagasituleku aeg T(x) ja Poincaré kujutus P(x).
    return [t[first], xsol[first]]

Selle funktsiooni abil võime nüüd $T(x)$ ja $P(x)$ arvutada.

In [36]:
xx = np.linspace(0.01, np.pi - 0.01, 500)
xx = np.concatenate((np.flipud(-xx), xx), 0)
res = np.array([poincd(x) for x in xx])
tt = res[:, 0]
pp = res[:, 1]

Joonistame kõigepealt $T(x)$. Joonis näitab, et tipu $x = \pm\pi$ ümbruses tagasitulekuaeg hajub, $T(x) \to \infty$. See on sellepärast, et sadula $(\pm\pi, 0)$ ümbruses liigub pendel aeglaselt, ja võnkumine võtab kauem aega.

In [37]:
fig = plt.figure()
ax = fig.add_subplot(111)
ax.set_xlabel('x')
ax.set_ylabel('T(x)')
ax.plot(xx, tt)
ax.grid()
plt.show()

Järgmisena joonistame $P(x)$. Paneme tähele, et kehtib $0 < |P(x)| < |x|$ ja $\mathrm{sgn} P(x) = -\mathrm{sgn} x$, kuna pendel kaotab energiat, ja Poincaré lõige on valitud sellisena, et pendli asukoht $x$ vahetab märki igal võnkumisel.

In [38]:
fig = plt.figure()
ax = fig.add_subplot(111)
ax.set_xlabel('x')
ax.set_ylabel('P(x)')
ax.plot(xx, pp)
ax.grid()
plt.show()