Skip to content

Pendulum with a vertically or horizontally oscillating pivot

A pendulum whose suspension point oscillates vertically or horizontally, with acceleration ε g cos Ω t.

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

Select this system

from skcomplex.systems import create_simulation, create_system, describe

system = create_system('parametric_pendulum')
entry = describe('parametric_pendulum')
trajectory = create_simulation('parametric_pendulum').prepare().run().trajectory_

Factory

@register_system(
    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,
    direction: str = 'vertical',
) -> System:
    """Return a pendulum whose suspension point oscillates vertically or horizontally, with acceleration `ε g cos Ω t`.

    Vertically, the acceleration of the suspension point modulates gravity, `θ'' + c θ' + (g / l)(1 + ε cos Ω t) sin θ
    = 0`. Small oscillations grow when `Ω` is close to twice the natural frequency `sqrt(g / l)`, the parametric
    resonance of a swing, while a pendulum at rest stays at rest. Horizontally, the acceleration pushes the bob
    sideways, `θ'' + c θ' + (g / l)(sin θ + ε cos Ω t cos θ) = 0`, which for small angles is a harmonic oscillator under
    an external force. Its amplitude grows linearly in time, as `ε Ω t / 2`, when `Ω` equals the natural frequency,
    the linear resonance, and the pendulum starts to swing even from rest.

    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.

        direction:
            The direction in which the suspension point oscillates, `'vertical'` or `'horizontal'`.

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

    Examples:
        >>> from skcomplex.systems.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', direction=direction)],
    )

Components

Elements:

  • ParametricPendulumElement, from Element

Interactions:

  • ParametricPendulumInteraction, from EulerLagrangeInteraction
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.systems.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 pendulum whose suspension point oscillates vertically or horizontally, an Euler-Lagrange equation.

    The suspension point accelerates by `ε g cos Ω t`. Vertically, the Lagrangian per moment of inertia is
    `θ'² / 2 + (g / l)(1 + ε cos Ω t) cos θ`, and the acceleration modulates gravity, a parametric forcing.
    Horizontally, it is `θ'² / 2 + (g / l)(cos θ - ε cos Ω t sin θ)`, and the acceleration pushes the bob sideways, an
    external forcing. The Rayleigh dissipation `c θ'² / 2` adds the friction, so the angle follows
    `θ'' + c θ' + (g / l)(1 + ε cos Ω t) sin θ = 0` or `θ'' + c θ' + (g / l)(sin θ + ε cos Ω t cos θ) = 0`, with the
    parameters of the target.

    Args:
        label:
            The label of the interaction.

        target:
            The path of the pendulum.

        direction:
            The direction in which the suspension point oscillates, `'vertical'` or `'horizontal'`.

    Raises:
        ValidationError:
            The direction is unknown.

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

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

    def _check(self) -> None:
        """Raise an error for an unknown direction of the suspension point.

        Raises:
            ValidationError:
                The direction is unknown.
        """
        super()._check()
        if self.direction not in DIRECTIONS:
            error_msg = f'Unknown direction {self.direction!r}. Use one of {", ".join(DIRECTIONS)}.'
            raise ValidationError(error_msg)

    @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
        drive = p.epsilon * cos(p.Omega * namespace.t)
        if self.direction == 'horizontal':
            return derivative(p.theta, 1) ** 2 / 2 + p.g / p.l * (cos(p.theta) - drive * sin(p.theta))
        return derivative(p.theta, 1) ** 2 / 2 + p.g / p.l * (1 + drive) * 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)