Parametric pendulum

A pendulum whose suspension point oscillates vertically, θ'' + c θ' + (g / l)(1 + ε cos Ω t) sin θ = 0.

• Name: parametric_pendulum
• Field: physics/mechanics
• Formalism: ordinary differential equations
• References: Dilão 2023

Factory

@register_model(
    field='physics/mechanics',
    tags=('mechanics', 'ode', 'second order', 'forcing', 'parametric resonance'),
    references=('Dilão 2023',),
    defaults={'until': 30.0, 'dt': 0.02},
)
def parametric_pendulum(
    theta: float = 0.1,
    omega: float = 0.0,
    *,
    g: float = 9.81,
    l: float = 1.0,
    epsilon: float = 0.2,
    Omega: float = 2 * math.sqrt(9.81),
    c: float = 0.0,
) -> System:
    """Return a pendulum whose suspension point oscillates vertically, `θ'' + c θ' + (g / l)(1 + ε cos Ω t) sin θ = 0`.

    The vertical acceleration of the suspension point modulates gravity. Small oscillations grow when `Ω` is close to
    twice the natural frequency `sqrt(g / l)`, the parametric resonance of a swing.

    Args:
        theta:
            The initial angle.

        omega:
            The initial angular velocity.

        g:
            The gravitational acceleration.

        l:
            The length of the rod.

        epsilon:
            The amplitude of the acceleration of the suspension point, relative to `g`.

        Omega:
            The angular frequency of the suspension point.

        c:
            The friction coefficient.

    Returns:
        system:
            The model, with the angle `swing.theta` of order 2.

    Examples:
        >>> from skcomplex.models.physics.mechanics import parametric_pendulum
        >>> parametric_pendulum(epsilon=0.5).elements[0].epsilon
        0.5
    """
    return System(
        elements=[
            ParametricPendulumElement(
                'swing',
                theta=theta,
                theta_d1=omega,
                g=g,
                l=l,
                epsilon=epsilon,
                Omega=Omega,
                c=c,
            ),
        ],
        interactions=[ParametricPendulumInteraction('motion', target='swing')],
    )

Components

class ParametricPendulumElement(Element):
    """Pendulum whose suspension point oscillates vertically, with the angle `theta` of order 2.

    The motion comes from `ParametricPendulumInteraction`.

    Args:
        label:
            The label of the pendulum.

        visual:
            How the pendulum is drawn, or `None` for none.

        theta:
            The initial angle.

        theta_d1:
            The initial angular velocity.

        g:
            The gravitational acceleration.

        l:
            The length of the rod.

        epsilon:
            The amplitude of the acceleration of the suspension point, relative to `g`.

        Omega:
            The angular frequency of the suspension point.

        c:
            The friction coefficient.

    Examples:
        >>> from skcomplex.models.physics.mechanics import ParametricPendulumElement
        >>> ParametricPendulumElement('swing', epsilon=0.5).epsilon
        0.5
    """

    def __init__(
        self,
        label: str,
        visual: VisualSpec | None = None,
        *,
        theta: Annotated[float, StateSpec(period=2 * pi, order=2)] = 0.1,
        theta_d1: float = 0.0,
        g: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 9.81,
        l: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 1.0,
        epsilon: Annotated[float, ParameterSpec()] = 0.2,
        Omega: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 2 * math.sqrt(9.81),
        c: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 0.0,
    ) -> None:
        super().__init__(label, visual=visual)
        self.theta = theta
        self.theta_d1 = theta_d1
        self.g = g
        self.l = l
        self.epsilon = epsilon
        self.Omega = Omega
        self.c = c

class ParametricPendulumInteraction(EulerLagrangeInteraction):
    """Motion of a parametric pendulum, the Euler-Lagrange equation of `θ'² / 2 + (g / l)(1 + ε cos Ω t) cos θ`.

    The Rayleigh dissipation `c θ'² / 2` adds the friction, so the angle follows
    `θ'' + c θ' + (g / l)(1 + ε cos Ω t) sin θ = 0`, with the parameters of the target.

    Args:
        label:
            The label of the interaction.

        target:
            The path of the pendulum.

    Examples:
        >>> from skcomplex.models.physics.mechanics import ParametricPendulumInteraction
        >>> ParametricPendulumInteraction('swing', target='pendulum').states
        'theta'
    """

    def __init__(self, label: str, *, target: Annotated[str, RoleSpec()]) -> None:
        super().__init__(label, target=target, states='theta', dissipation=self._get_friction)

    @staticmethod
    def _get_friction(namespace: SimpleNamespace) -> Expression:
        """Return the Rayleigh dissipation `c θ'² / 2` of the friction."""
        return namespace.target.c * derivative(namespace.target.theta, 1) ** 2 / 2

    def _get_lagrangian(self, namespace: SimpleNamespace) -> Expression:
        """Return the Lagrangian of the parametric pendulum per moment of inertia.

        Args:
            namespace:
                The namespace of the symbols the Lagrangian can use.

        Returns:
            lagrangian:
                The Lagrangian.
        """
        p = namespace.target
        gravity = p.g / p.l * (1 + p.epsilon * cos(p.Omega * namespace.t))
        return derivative(p.theta, 1) ** 2 / 2 + gravity * cos(p.theta)

Equations

swing.theta: diff(theta_swing(t), t, t) == -(epsilon_swing*cos(Omega_swing*t) + 1)*g_swing*sin(theta_swing(t))/l_swing - c_swing*diff(theta_swing(t), t)