Skip to content

Double pendulum

A double pendulum, a pendulum hanging from the bob of another, written as its two Lagrange equations.

  • Name: double_pendulum
  • Field: physics/mechanics
  • Formalism: ordinary differential equations
  • References: Ryzhov et al. 2021; Argun et al. 2021

Factory

@register_model(
    field='physics/mechanics',
    tags=('mechanics', 'ode', 'second order', 'chaos', 'lagrangian'),
    references=('Ryzhov et al. 2021', 'Argun et al. 2021'),
    defaults={'until': 20.0, 'dt': 0.02},
)
def double_pendulum(
    theta1: float = math.pi / 2,
    theta2: float = math.pi / 2,
    omega1: float = 0.0,
    omega2: float = 0.0,
    *,
    m1: float = 1.0,
    m2: float = 1.0,
    l1: float = 1.0,
    l2: float = 1.0,
    g: float = 9.81,
) -> System:
    """Return a double pendulum, a pendulum hanging from the bob of another, written as its two Lagrange equations.

    The equations `(m1 + m2) l1 θ1'' + m2 l2 θ2'' cos(θ1 - θ2) + m2 l2 θ2'² sin(θ1 - θ2) + (m1 + m2) g sin θ1 = 0`
    and `l2 θ2'' + l1 θ1'' cos(θ1 - θ2) - l1 θ1'² sin(θ1 - θ2) + g sin θ2 = 0` are solved jointly for the angular
    accelerations. The motion is chaotic at large amplitudes.

    Args:
        theta1:
            The initial angle of the upper rod from the vertical, in radians.

        theta2:
            The initial angle of the lower rod from the vertical, in radians.

        omega1:
            The initial angular velocity of the upper rod.

        omega2:
            The initial angular velocity of the lower rod.

        m1:
            The mass of the upper bob.

        m2:
            The mass of the lower bob.

        l1:
            The length of the upper rod.

        l2:
            The length of the lower rod.

        g:
            The gravitational acceleration.

    Returns:
        system:
            The model, with the angles `pendulum.theta1` and `pendulum.theta2` of order 2.

    Examples:
        >>> from skcomplex.models.physics.mechanics import double_pendulum
        >>> double_pendulum(m2=2.0).elements[0].m2
        2.0
    """
    return System(
        elements=[
            DoublePendulumElement(
                'pendulum',
                theta1=theta1,
                theta2=theta2,
                theta1_d1=omega1,
                theta2_d1=omega2,
                m1=m1,
                m2=m2,
                l1=l1,
                l2=l2,
                g=g,
            ),
        ],
        interactions=[DoublePendulumInteraction('motion', target='pendulum')],
    )

Components

class DoublePendulumElement(Element):
    """Double pendulum, a bob of mass `m2` on a rod of length `l2` hanging from a bob of mass `m1` on a rod `l1`.

    The angles `theta1` and `theta2` of the rods from the vertical have order 2. The motion comes from
    `DoublePendulumInteraction`.

    Args:
        label:
            The label of the double pendulum.

        visual:
            How the double pendulum is drawn, or `None` for its class visual.

        theta1:
            The initial angle of the upper rod.

        theta1_d1:
            The initial angular velocity of the upper rod.

        theta2:
            The initial angle of the lower rod.

        theta2_d1:
            The initial angular velocity of the lower rod.

        m1:
            The mass of the upper bob.

        m2:
            The mass of the lower bob.

        l1:
            The length of the upper rod.

        l2:
            The length of the lower rod.

        g:
            The gravitational acceleration.

    Examples:
        >>> from skcomplex.models.physics.mechanics import DoublePendulumElement
        >>> DoublePendulumElement('pendulum', m2=2.0).m2
        2.0
    """

    visual = VisualSpec('bodies', {'positions': _double_bobs, 'pivot': (0.0, 0.0), 'rods': True, 'traced': [1]})

    def __init__(
        self,
        label: str,
        visual: VisualSpec | None = None,
        *,
        theta1: Annotated[float, StateSpec(period=2 * pi, order=2)] = math.pi / 2,
        theta1_d1: float = 0.0,
        theta2: Annotated[float, StateSpec(period=2 * pi, order=2)] = math.pi / 2,
        theta2_d1: float = 0.0,
        m1: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 1.0,
        m2: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 1.0,
        l1: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 1.0,
        l2: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 1.0,
        g: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 9.81,
    ) -> None:
        super().__init__(label, visual=visual)
        self.theta1 = theta1
        self.theta1_d1 = theta1_d1
        self.theta2 = theta2
        self.theta2_d1 = theta2_d1
        self.m1 = m1
        self.m2 = m2
        self.l1 = l1
        self.l2 = l2
        self.g = g

class DoublePendulumInteraction(EulerLagrangeInteraction):
    """Motion of a double pendulum, the Euler-Lagrange equations of its kinetic energy and gravity.

    The Lagrangian is `(m1 + m2) l1² θ1'² / 2 + m2 l2² θ2'² / 2 + m2 l1 l2 θ1' θ2' cos(θ1 - θ2) + (m1 + m2) g l1 cos θ1
    + m2 g l2 cos θ2`, with the masses, lengths, and gravity of the target.

    Args:
        label:
            The label of the interaction.

        target:
            The path of the double pendulum.

    Examples:
        >>> from skcomplex.models.physics.mechanics import DoublePendulumInteraction
        >>> DoublePendulumInteraction('motion', target='pendulum').states
        ('theta1', 'theta2')
    """

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

    def _get_lagrangian(self, namespace: SimpleNamespace) -> Expression:
        """Return the Lagrangian of the double pendulum.

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

        Returns:
            lagrangian:
                The Lagrangian.
        """
        p = namespace.target
        w1, w2 = derivative(p.theta1, 1), derivative(p.theta2, 1)
        kinetic = (p.m1 + p.m2) * p.l1**2 * w1**2 / 2 + p.m2 * p.l2**2 * w2**2 / 2
        kinetic += p.m2 * p.l1 * p.l2 * w1 * w2 * cos(p.theta1 - p.theta2)
        return kinetic + (p.m1 + p.m2) * p.g * p.l1 * cos(p.theta1) + p.m2 * p.g * p.l2 * cos(p.theta2)

Equations

pendulum.theta1: diff(theta1_pendulum(t), t, t) == (l1_pendulum*l2_pendulum*m2_pendulum*sin(theta1_pendulum(t) - theta2_pendulum(t))*diff(theta1_pendulum(t), t)^2 - g_pendulum*l2_pendulum*m2_pendulum*sin(theta2_pendulum(t)) + (l1_pendulum*l2_pendulum*m2_pendulum*sin(theta1_pendulum(t) - theta2_pendulum(t))*diff(theta2_pendulum(t), t)^2 + g_pendulum*l1_pendulum*(m1_pendulum + m2_pendulum)*sin(theta1_pendulum(t)))*l2_pendulum*m2_pendulum*cos(theta1_pendulum(t) - theta2_pendulum(t))/(l1_pendulum*(m1_pendulum + m2_pendulum)))*l2_pendulum*m2_pendulum*cos(theta1_pendulum(t) - theta2_pendulum(t))/((l2_pendulum^2*m2_pendulum^2*cos(theta1_pendulum(t) - theta2_pendulum(t))^2/(m1_pendulum + m2_pendulum) - l2_pendulum^2*m2_pendulum)*l1_pendulum*(m1_pendulum + m2_pendulum)) - (l1_pendulum*l2_pendulum*m2_pendulum*sin(theta1_pendulum(t) - theta2_pendulum(t))*diff(theta2_pendulum(t), t)^2 + g_pendulum*l1_pendulum*(m1_pendulum + m2_pendulum)*sin(theta1_pendulum(t)))/(l1_pendulum^2*(m1_pendulum + m2_pendulum))
pendulum.theta2: diff(theta2_pendulum(t), t, t) == -(l1_pendulum*l2_pendulum*m2_pendulum*sin(theta1_pendulum(t) - theta2_pendulum(t))*diff(theta1_pendulum(t), t)^2 - g_pendulum*l2_pendulum*m2_pendulum*sin(theta2_pendulum(t)) + (l1_pendulum*l2_pendulum*m2_pendulum*sin(theta1_pendulum(t) - theta2_pendulum(t))*diff(theta2_pendulum(t), t)^2 + g_pendulum*l1_pendulum*(m1_pendulum + m2_pendulum)*sin(theta1_pendulum(t)))*l2_pendulum*m2_pendulum*cos(theta1_pendulum(t) - theta2_pendulum(t))/(l1_pendulum*(m1_pendulum + m2_pendulum)))/(l2_pendulum^2*m2_pendulum^2*cos(theta1_pendulum(t) - theta2_pendulum(t))^2/(m1_pendulum + m2_pendulum) - l2_pendulum^2*m2_pendulum)