%%html
<style>
.output_wrapper, .output {
    height:auto !important;
    max-height:10000px;
}
.output_scroll {
    box-shadow:none !important;
    webkit-box-shadow:none !important;
}
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>

Mittelineaarne pendel ja sõltuvus algtingimustest

See tööleht näitab, kuidas saab simuleerida matemaatilist pendli koos välisjõuga, mis tekitab mittelineaarset võnkumist. Pendli liikumisvõrrandite lahendamiseks kasutame SciPy ja NumPy, tulemuste joonistamiseks kasutame PyPlot. Selleks et näha kuidas pendli liikumine sõltib algtingimustest, valime mitu algtingimust ja joonistame lahandeid ühises joonises.

%matplotlib inline

from scipy.integrate import odeint
import matplotlib.pyplot as plt
import numpy as np
from mpl_toolkits.mplot3d import Axes3D
import multiprocessing as mp
from itertools import product

Pendli definitsioon

Kui matemaatilise pendli pikkus on $l$, mass on $m$ ja ta liigub gravitatsiooniväljas mille vabalangemiskiirendus on $g$, siis kirjeldab tema liikumist teist järku diferentsiaalvõrrand

$$ml\ddot{\theta} = -mg\sin\theta\,,$$

kus $\theta$ on nurk pendli asukoha ja alumise positsiooni vahel. Kui lisaks tekib ka hõõrdejõud, mis sõltub lineaarselt pendli kiirusest ja mille hõõrdetegur on $h$, siis tuleb liikumisvõrrandisse vastav liige,

$$ml\ddot{\theta} = -mg\sin\theta - hl\dot{\theta}\,.$$

Kui lisaks sellele rakendub ka perioodiline väline jõud, mille nurksagedus on $\Omega$ ja amplituud on $F$, tuleb veel üks liige juurde,

$$ml\ddot{\theta} = -mg\sin\theta - hl\dot{\theta} - F\cos(\Omega t)\,.$$

Simulatsiooni jaoks võib liikumisvõrrandit kõigepealt lihtsustada. Selleks defineerime

$$K = \frac{g}{l}\,, \quad H = \frac{h}{m}\,, \quad A = \frac{F}{ml}\,,$$

ja kirjutame liikumisvõrrandit kujul

$$\ddot{\theta} = -K\sin\theta - H\dot{\theta} - A\cos(\Omega t)\,.$$

See on teist järku diferentsiaalvõrrand. Selleks et kasutada dünaamiliste süsteemide teeoriat, on vaja tuletada esimest järku diferentsiaalvõrrandit. Selleks defineerime $\omega = \dot{\theta}$ ja saame süsteemi

\begin{align} \dot{\theta} &= \omega\,,\\ \dot{\omega} &= -K\sin\theta - H\omega - A\cos(\Omega t)\,. \end{align}

Kuna dünaamika sõltub ka ajast $t$, on süsteem mitteautonoomne. Autonoomset süsteemi saame näiteks definitsiooni $\varphi = \Omega t$ abil, millest järeldub

\begin{align} \dot{\theta} &= \omega\,,\\ \dot{\omega} &= -K\sin\theta - H\omega - A\cos(\varphi)\,,\\ \dot{\varphi} &= \Omega\,, \end{align}

koos algtingimusega $\varphi(0) = 0$.

Pendli numbriline simulatsioon

Simulatsiooni jaoks on siiski lihtsam kasutada mitteautonoomset süsteemi. Selleks defineerime kõigepealt funktsiooni, mis kirjeldab süsteemi dünaamikat:

def pend(x, t, kk, hh, aa, ww):
    theta, omega = x
    dxdt = [omega, -kk * np.sin(theta) - hh * omega - aa * np.cos(ww * t)]
    return dxdt

Muutuja x kirjeldab pendli olekut $(\theta, \omega)$. Võrrandite lahendamise jaoks on meil lisaks ka algtingimusi vaja. Kuna me soovime võrrandeid lahendada korraga mitme algtingimuse jaoks, defineerime algtingimuste keskpunkti x0 ja nende vahelist kaugust dx0. Muutuja count määrab mitu algtingimuste väärtust me kasutame iga dimensiooni kohta.

x0 = [0.0, 0.0]
dx0 = [0.1, 0.1]
count = 8

theta0 = np.linspace(x0[0] - (count - 1) * dx0[0] / 2, x0[0] + (count - 1) * dx0[0] / 2, count)
omega0 = np.linspace(x0[1] - (count - 1) * dx0[1] / 2, x0[1] + (count - 1) * dx0[1] / 2, count)

ini = list(product(theta0, omega0))

Järgmisena defineerime simulatsiooni aega. time on simulatsiooni kestus $t_{\text{max}}$, ehk me lahendame võrrandeid ajavahemikus $t \in [0, t_{\text{max}}]$. rate on simulatsioonisammude arv ajaühiku kohta. Me kasutame NumPy funktsiooni linspace ajaväärtuste vektori loomiseks.

time = 25
rate = 100

samples = time * rate + 1
t = np.linspace(0, time, samples)

Lõpuks defineerime funktsiooni, mis lahendab liikumisvõrrandeid. Selleks kasutame SciPy funktsiooni odeint. Tulemusi näitame mitmes joonises, kus on joonistatud erinevad muutujad. Paralleelarvutuse rakendamiseks kasutame objektklassi.

class pendint(object):
    def __init__(self, time, kk, hh, aa, ww):
        self.time = time
        self.param = (kk, hh, aa, ww)
        
    def __call__(self, y0):
        return odeint(pend, y0, self.time, args = self.param)

def sim(kk, hh, aa, ww):
    # Lahendame liikumisvõrrandeid.
    intobj = pendint(t, kk, hh, aa, ww)
    with mp.Pool() as pool:
        sol = np.array(pool.map(intobj, ini)).reshape(count, count, samples, 2)
    
    # Joonis: aeg ja asukoht
    fig = plt.figure(figsize=(16, 12))
    ax = fig.add_subplot(111)
    for s1 in sol:
        for s0 in s1:
            ax.plot(t, s0[:, 0])
    ax.set_xlabel("t")
    ax.set_ylabel("$\\theta$")
    ax.grid()
    plt.show()
    
    # Joonis: aeg ja kiirus
    fig = plt.figure(figsize=(16, 12))
    ax = fig.add_subplot(111)
    for s1 in sol:
        for s0 in s1:
            ax.plot(t, s0[:, 1])
    ax.set_xlabel("t")
    ax.set_ylabel("$\\omega$")
    ax.grid()
    plt.show()
    
    # Joonis: asukoht ja kiirus
    fig = plt.figure(figsize=(16, 12))
    ax = fig.add_subplot(111)
    for s1 in sol:
        for s0 in s1:
            ax.plot(s0[:, 0], s0[:, 1])
    ax.set_xlabel("$\\theta$")
    ax.set_ylabel("$\\omega$")
    ax.grid()
    plt.show()
    
    # Joonis: ruumiline joonis
    fig = plt.figure(figsize=(16, 12))
    ax = fig.add_subplot(111, projection = '3d')
    for s1 in sol:
        for s0 in s1:
            ax.plot(t, s0[:, 0], s0[:, 1])
    ax.set_xlabel("t")
    ax.set_ylabel("$\\theta$")
    ax.set_zlabel("$\\omega$")
    plt.show()

Simulatsiooni tulemused

Hõõrdejõuga pendel ilma välisjõuta

Kui välisjõudu ei ole, summutab hõõrdejõud pendli võnkumist. Trajektoorid koonduvad stabiilsele püsipunktile $\theta = 0, \omega = 0$.

sim(5.0, 1.0, 0.0, 0.0)

sim(5.0, 0.0, 0.0, 0.0)

sim(5.0, 1.0, 1.0, 1.2)

sim(5.0, 1.0, 5.0, 1.2)

sim(5.0, 1.0, 10.0, 1.2)