Skip to content

Pendulum¤

Pendulum ¤

Generate time series data for a pendulum.

We describe a generic pendulum system by the Lagrangian action $$ S_L[\theta] = I \int_{t_0}^{t_1} \mathbb{d}t \left\{\frac{1}{2} \dot\theta^2 + \omega_0^2 \cos\theta \right\}\,, $$ where \(\theta\) is the angle from the vertical to the pendulum; \(I\) is the inertia parameter introduced for dimensional reasons, and \(\omega_0\) the frequency parameter.

Details are collected in the tutorial.

Source code in hamilflow/models/pendulum.py 
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
class Pendulum:
    r"""Generate time series data for a pendulum.

    We describe a generic pendulum system by the Lagrangian action
    $$
    S_L\[\theta\] = I \int_{t_0}^{t_1} \mathbb{d}t
    \left\\{\frac{1}{2} \dot\theta^2 + \omega_0^2 \cos\theta \right\\}\,,
    $$
    where $\theta$ is the _angle_ from the vertical to the pendulum;
    $I$ is the _inertia parameter_ introduced for dimensional reasons,
    and $\omega_0$ the _frequency parameter_.

    Details are collected in the tutorial.
    """

    def __init__(
        self,
        system: int | float | dict[str, int | float],
        initial_condition: int | float | dict[str, int | float],
    ) -> None:
        if isinstance(system, (float, int)):
            system = {"omega0": system}
        if isinstance(initial_condition, (float, int)):
            initial_condition = {"theta0": initial_condition}
        self.system = PendulumSystem.model_validate(system)
        self.initial_condition = PendulumIC.model_validate(initial_condition)

    @cached_property
    def definition(self) -> dict[str, float]:
        """Model params and initial conditions defined as a dictionary."""
        return dict(
            system=self.system.model_dump(),
            initial_condition=self.initial_condition.model_dump(),
        )

    @property
    def omega0(self) -> float:
        return self.system.omega0

    @property
    def _k(self) -> float:
        return self.initial_condition.k

    @property
    def _math_m(self) -> float:
        return self._k**2

    @cached_property
    def freq(self) -> float:
        r"""Frequency.

        :return: $\frac{\pi}{2K(k^2)}\omega_0$, where
        $K(m)$ is [Legendre's complete elliptic integral of the first kind](https://dlmf.nist.gov/19.2#E8)
        """
        return math.pi * self.omega0 / (2 * ellipk(self._math_m))

    @cached_property
    def period(self) -> float:
        r"""Period.

        :return: $\frac{4K(k^2)}{\omega_0}$, where
        $K(m)$ is [Legendre's complete elliptic integral of the first kind](https://dlmf.nist.gov/19.2#E8)
        """
        return 4 * ellipk(self._math_m) / self.omega0

    def _math_u(self, t: ArrayLike) -> np.ndarray[float]:
        return self.omega0 * np.asarray(t)

    def u(self, t: ArrayLike) -> np.ndarray[float]:
        r"""The convenient generalised coordinate $u$,
        $\sin u \coloneqq \frac{\sin\frac{\theta}{2}}{\sin\frac{\theta_0}{2}}$.

        :param t: time
        :return: $u(t) = \mathrm{am}\!\big(\omega_0 t + K(k^2), k^2\big)$, where
        $\mathrm{am}(x, k)$ is [Jacobi's amplitude function](https://dlmf.nist.gov/22.16#E1),
        $K(m)$ is [Legendre's complete elliptic integral of the first kind](https://dlmf.nist.gov/19.2#E8)
        """
        _, _, _, ph = ellipj(self._math_u(t) + ellipk(self._math_m), self._math_m)

        return ph

    def theta(self, t: ArrayLike) -> np.ndarray[float]:
        r"""Angle $\theta$.

        :param t: time
        :return: $\theta(t) = 2\arcsin\!\big(k\cdot\mathrm{cd}(\omega_0 t, k^2)\big)$, where
        $\mathrm{cd}(z, k)$ is a [Jacobian elliptic function](https://dlmf.nist.gov/22.2#E8)
        """
        _, cn, dn, _ = ellipj(self._math_u(t), self._math_m)

        return 2 * np.arcsin(cn / dn * self._k)

    def __call__(self, n_periods: int, n_samples_per_period: int) -> pd.DataFrame:
        time_delta = self.period / n_samples_per_period
        time_steps = np.arange(0, n_periods * n_samples_per_period) * time_delta

        thetas = self.theta(time_steps)
        us = self.u(time_steps)

        return pd.DataFrame(dict(t=time_steps, x=thetas, u=us))

definition: dict[str, float] cached property ¤

Model params and initial conditions defined as a dictionary.

freq: float cached property ¤

Frequency.

Returns:

Type Description
float

\(\frac{\pi}{2K(k^2)}\omega_0\), where \(K(m)\) is Legendre's complete elliptic integral of the first kind

period: float cached property ¤

Period.

Returns:

Type Description
float

\(\frac{4K(k^2)}{\omega_0}\), where \(K(m)\) is Legendre's complete elliptic integral of the first kind

theta(t) ¤

Angle \(\theta\).

Parameters:

Name Type Description Default
t ArrayLike

time

 required

Returns:

Type Description
ndarray[float]

\(\theta(t) = 2\arcsin\!\big(k\cdot\mathrm{cd}(\omega_0 t, k^2)\big)\), where \(\mathrm{cd}(z, k)\) is a Jacobian elliptic function
Source code in hamilflow/models/pendulum.py 
121
122
123
124
125
126
127
128
129
130
def theta(self, t: ArrayLike) -> np.ndarray[float]:
    r"""Angle $\theta$.

    :param t: time
    :return: $\theta(t) = 2\arcsin\!\big(k\cdot\mathrm{cd}(\omega_0 t, k^2)\big)$, where
    $\mathrm{cd}(z, k)$ is a [Jacobian elliptic function](https://dlmf.nist.gov/22.2#E8)
    """
    _, cn, dn, _ = ellipj(self._math_u(t), self._math_m)

    return 2 * np.arcsin(cn / dn * self._k)

u(t) ¤

The convenient generalised coordinate \(u\), \(\sin u \coloneqq \frac{\sin\frac{\theta}{2}}{\sin\frac{\theta_0}{2}}\).

Parameters:

Name Type Description Default
t ArrayLike

time

 required

Returns:

Type Description
ndarray[float]

\(u(t) = \mathrm{am}\!\big(\omega_0 t + K(k^2), k^2\big)\), where \(\mathrm{am}(x, k)\) is Jacobi's amplitude function, \(K(m)\) is Legendre's complete elliptic integral of the first kind
Source code in hamilflow/models/pendulum.py 
108
109
110
111
112
113
114
115
116
117
118
119
def u(self, t: ArrayLike) -> np.ndarray[float]:
    r"""The convenient generalised coordinate $u$,
    $\sin u \coloneqq \frac{\sin\frac{\theta}{2}}{\sin\frac{\theta_0}{2}}$.

    :param t: time
    :return: $u(t) = \mathrm{am}\!\big(\omega_0 t + K(k^2), k^2\big)$, where
    $\mathrm{am}(x, k)$ is [Jacobi's amplitude function](https://dlmf.nist.gov/22.16#E1),
    $K(m)$ is [Legendre's complete elliptic integral of the first kind](https://dlmf.nist.gov/19.2#E8)
    """
    _, _, _, ph = ellipj(self._math_u(t) + ellipk(self._math_m), self._math_m)

    return ph

PendulumIC ¤

Bases: BaseModel

The initial condition for a pendulum.

Parameters:

Name Type Description Default
theta0

\(-\frac{\pi}{2} \le \theta_0 \le \frac{\pi}{2}\), the initial angle

 required
Source code in hamilflow/models/pendulum.py 
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
class PendulumIC(BaseModel):
    r"""The initial condition for a pendulum.

    :param theta0: $-\frac{\pi}{2} \le \theta_0 \le \frac{\pi}{2}$, the
    initial angle
    """

    theta0: float = Field(ge=-math.pi / 2, le=math.pi / 2, frozen=True)

    @computed_field  # type: ignore[misc]
    @cached_property
    def k(self) -> float:
        r"""A convenient number for elliptic functions.

        :return: $\sin\frac{\theta_0}{2}$
        """
        return math.sin(self.theta0 / 2)

k: float cached property ¤

A convenient number for elliptic functions.

Returns:

Type Description
float

\(\sin\frac{\theta_0}{2}\)

PendulumSystem ¤

Bases: BaseModel

The params for the pendulum.

Parameters:

Name Type Description Default
omega0

\(\omega_0 \coloneqq \sqrt{\frac{U}{I}} > 0\), frequency parameter

 required
Source code in hamilflow/models/pendulum.py 
11
12
13
14
15
16
17
18
class PendulumSystem(BaseModel):
    r"""The params for the pendulum.

    :param omega0: $\omega_0 \coloneqq \sqrt{\frac{U}{I}} > 0$, frequency
    parameter
    """

    omega0: float = Field(gt=0, frozen=True)