Skip to content

Huygens clocks and pendulums on coupled supports

Pendulum clocks on supports joined by springs, which synchronize through the motion of the supports.

  • Name: huygens_clocks
  • Field: physics/mechanics
  • Formalism: ordinary differential equations
  • References: Huygens 1673; Dilão 2023; Dilão 2023

Factory

@register_model(
    field='physics/mechanics',
    tags=('mechanics', 'ode', 'second order', 'synchronization', 'limit cycle'),
    references=('Huygens 1673', 'Dilão 2023', 'Dilão 2023'),
    defaults={'until': 20.0, 'dt': 0.05},
)
def huygens_clocks(
    n_clocks: int = 2,
    *,
    theta: float = 0.2,
    k: float = 1.0,
    m: float = 1.0,
    l: float = 1.0,
    M: float = 1.0,
    g: float = 9.8,
    rho: float = 0.1,
    lam: float = 0.1,
    theta_c: float = 0.1,
) -> System:
    """Return pendulum clocks on supports joined by springs, which synchronize through the motion of the supports.

    Every pendulum hangs from a support of mass `M` that moves horizontally with damping `rho`, and neighbouring
    supports are joined by springs of stiffness `k`. A pendulum follows `m l θ'' + f(θ) θ' + m g sin θ = -m x'' cos θ`
    and its support `(M + m) x'' + 2 rho x' + m l θ'' cos θ = m l θ'² sin θ + spring forces`. The escapement
    `f(θ) = 2 λ sgn(|θ| - θ_c)` is `-2 λ` for `|θ| < θ_c`, pumping energy in, and `2 λ` beyond, so every clock has a
    limit cycle. With `λ = 0` the clocks are plain pendulums. Two identical clocks synchronize in antiphase. The
    first clock starts at `theta`, away from the switch of the escapement at `θ_c`, and the others at rest.

    Args:
        n_clocks:
            The number of clocks in the row.

        theta:
            The initial angle of the first clock. The others start at rest.

        k:
            The stiffness of the springs between neighbouring supports.

        m:
            The mass of every bob.

        l:
            The length of every pendulum.

        M:
            The mass of every support.

        g:
            The gravitational acceleration.

        rho:
            The damping of the supports.

        lam:
            The strength of the escapement.

        theta_c:
            The angle below which the escapement pumps energy in.

    Returns:
        system:
            The model, with the angles `clock_<i>.theta` and the support positions `clock_<i>.x`, all of order 2.

    Examples:
        >>> from skcomplex.models.physics.mechanics import huygens_clocks
        >>> len(huygens_clocks(3).interactions)
        5
    """
    clocks = [
        ClockPendulumElement(
            f'clock_{index}',
            theta=theta if index == 0 else 0.0,
            m=m,
            l=l,
            M=M,
            g=g,
            rho=rho,
            lam=lam,
            theta_c=theta_c,
        )
        for index in range(n_clocks)
    ]
    springs = [
        SupportSpringInteraction(f'spring_{index}', first=f'clock_{index}', second=f'clock_{index + 1}', k=k)
        for index in range(n_clocks - 1)
    ]
    motions = [ClockPendulumInteraction(f'motion_{index}', target=f'clock_{index}') for index in range(n_clocks)]
    return System(elements=clocks, interactions=[*springs, *motions])

Components

class ClockPendulumElement(Element):
    """Pendulum clock, a pendulum of angle `theta` hanging from a support at the horizontal position `x`.

    The angle and the position have order 2. The motion comes from `ClockPendulumInteraction`.

    Args:
        label:
            The label of the clock.

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

        theta:
            The initial angle of the pendulum.

        theta_d1:
            The initial angular velocity of the pendulum.

        x:
            The initial position of the support.

        x_d1:
            The initial velocity of the support.

        m:
            The mass of the bob.

        l:
            The length of the pendulum.

        M:
            The mass of the support.

        g:
            The gravitational acceleration.

        rho:
            The damping of the support.

        lam:
            The strength of the escapement.

        theta_c:
            The angle below which the escapement pumps energy in.

    Examples:
        >>> from skcomplex.models.physics.mechanics import ClockPendulumElement
        >>> ClockPendulumElement('clock', lam=0.2).lam
        0.2
    """

    def __init__(
        self,
        label: str,
        visual: VisualSpec | None = None,
        *,
        theta: Annotated[float, StateSpec(order=2)] = 0.2,
        theta_d1: float = 0.0,
        x: Annotated[float, StateSpec(order=2)] = 0.0,
        x_d1: float = 0.0,
        m: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 1.0,
        l: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 1.0,
        M: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 1.0,
        g: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 9.8,
        rho: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 0.1,
        lam: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 0.1,
        theta_c: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 0.1,
    ) -> None:
        super().__init__(label, visual=visual)
        self.theta = theta
        self.theta_d1 = theta_d1
        self.x = x
        self.x_d1 = x_d1
        self.m = m
        self.l = l
        self.M = M
        self.g = g
        self.rho = rho
        self.lam = lam
        self.theta_c = theta_c

class SupportSpringInteraction(EulerLagrangeInteraction):
    """Spring of stiffness `k` between the supports of two pendulum clocks.

    The spring adds the Lagrangian `-k (x2 - x1)² / 2` to the Euler-Lagrange equations of the supports, so every clock
    answers with the accelerations of its support and its pendulum through the mass matrix of its own Lagrangian.

    Args:
        label:
            The label of the interaction.

        first:
            The path of the first clock.

        second:
            The path of the second clock.

        k:
            The stiffness of the spring.

    Examples:
        >>> from skcomplex.models.physics.mechanics import SupportSpringInteraction
        >>> SupportSpringInteraction('spring', first='clock_0', second='clock_1', k=2.0).k
        2.0
    """

    def __init__(
        self,
        label: str,
        *,
        first: Annotated[str, RoleSpec()],
        second: Annotated[str, RoleSpec()],
        k: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 1.0,
    ) -> None:
        super().__init__(label, target=first, states='x')
        self.first = first
        self.second = second
        self.k = k

    def _get_lagrangian(self, namespace: SimpleNamespace) -> Expression:
        """Return minus the energy of the spring.

        Args:
            namespace:
                The namespace of the symbols the Lagrangian can use, with `first` and `second`.

        Returns:
            lagrangian:
                The Lagrangian of the spring.
        """
        return -namespace.k * (namespace.second.x - namespace.first.x) ** 2 / 2

class ClockPendulumInteraction(EulerLagrangeInteraction):
    """Motion of a pendulum clock on a moving support, the Euler-Lagrange equations with the escapement.

    The Lagrangian is `(M + m) x'² / 2 + m l x' θ' cos θ + m l² θ'² / 2 + m g l cos θ`. The Rayleigh dissipation
    `λ l sgn(|θ| - θ_c) θ'² + rho x'²` damps the support and gives the escapement `2 λ l sgn(|θ| - θ_c) θ'`, which pumps
    energy into the pendulum below the angle `θ_c` and takes it out beyond.

    Args:
        label:
            The label of the interaction.

        target:
            The path of the clock.

    Examples:
        >>> from skcomplex.models.physics.mechanics import ClockPendulumInteraction
        >>> ClockPendulumInteraction('motion', target='clock').states
        ('theta', 'x')
    """

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

    @staticmethod
    def _get_escapement(namespace: SimpleNamespace) -> Expression:
        """Return the Rayleigh dissipation of the escapement and of the damping of the support."""
        p = namespace.target
        escapement = p.lam * p.l * sign(abs(p.theta) - p.theta_c) * derivative(p.theta, 1) ** 2
        return escapement + p.rho * derivative(p.x, 1) ** 2

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

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

        Returns:
            lagrangian:
                The Lagrangian.
        """
        p = namespace.target
        v, w = derivative(p.x, 1), derivative(p.theta, 1)
        kinetic = (p.M + p.m) * v**2 / 2 + p.m * p.l * v * w * cos(p.theta) + p.m * p.l**2 * w**2 / 2
        return kinetic + p.m * p.g * p.l * cos(p.theta)

Equations

clock_0.theta: diff(theta_clock_0(t), t, t) == (g_clock_0*l_clock_0*m_clock_0*sin(theta_clock_0(t)) + 2*l_clock_0*lam_clock_0*sgn(-theta_c_clock_0 + abs(theta_clock_0(t)))*diff(theta_clock_0(t), t) + (l_clock_0*m_clock_0*sin(theta_clock_0(t))*diff(theta_clock_0(t), t)^2 - k_spring_0*(x_clock_0(t) - x_clock_1(t)) - 2*rho_clock_0*diff(x_clock_0(t), t))*l_clock_0*m_clock_0*cos(theta_clock_0(t))/(M_clock_0 + m_clock_0))/(l_clock_0^2*m_clock_0^2*cos(theta_clock_0(t))^2/(M_clock_0 + m_clock_0) - l_clock_0^2*m_clock_0)
clock_0.x: diff(x_clock_0(t), t, t) == -(g_clock_0*l_clock_0*m_clock_0*sin(theta_clock_0(t)) + 2*l_clock_0*lam_clock_0*sgn(-theta_c_clock_0 + abs(theta_clock_0(t)))*diff(theta_clock_0(t), t) + (l_clock_0*m_clock_0*sin(theta_clock_0(t))*diff(theta_clock_0(t), t)^2 - k_spring_0*(x_clock_0(t) - x_clock_1(t)) - 2*rho_clock_0*diff(x_clock_0(t), t))*l_clock_0*m_clock_0*cos(theta_clock_0(t))/(M_clock_0 + m_clock_0))*l_clock_0*m_clock_0*cos(theta_clock_0(t))/((l_clock_0^2*m_clock_0^2*cos(theta_clock_0(t))^2/(M_clock_0 + m_clock_0) - l_clock_0^2*m_clock_0)*(M_clock_0 + m_clock_0)) + (l_clock_0*m_clock_0*sin(theta_clock_0(t))*diff(theta_clock_0(t), t)^2 - k_spring_0*(x_clock_0(t) - x_clock_1(t)) - 2*rho_clock_0*diff(x_clock_0(t), t))/(M_clock_0 + m_clock_0)
clock_1.theta: diff(theta_clock_1(t), t, t) == (g_clock_1*l_clock_1*m_clock_1*sin(theta_clock_1(t)) + 2*l_clock_1*lam_clock_1*sgn(-theta_c_clock_1 + abs(theta_clock_1(t)))*diff(theta_clock_1(t), t) + (l_clock_1*m_clock_1*sin(theta_clock_1(t))*diff(theta_clock_1(t), t)^2 + k_spring_0*(x_clock_0(t) - x_clock_1(t)) - 2*rho_clock_1*diff(x_clock_1(t), t))*l_clock_1*m_clock_1*cos(theta_clock_1(t))/(M_clock_1 + m_clock_1))/(l_clock_1^2*m_clock_1^2*cos(theta_clock_1(t))^2/(M_clock_1 + m_clock_1) - l_clock_1^2*m_clock_1)
clock_1.x: diff(x_clock_1(t), t, t) == -(g_clock_1*l_clock_1*m_clock_1*sin(theta_clock_1(t)) + 2*l_clock_1*lam_clock_1*sgn(-theta_c_clock_1 + abs(theta_clock_1(t)))*diff(theta_clock_1(t), t) + (l_clock_1*m_clock_1*sin(theta_clock_1(t))*diff(theta_clock_1(t), t)^2 + k_spring_0*(x_clock_0(t) - x_clock_1(t)) - 2*rho_clock_1*diff(x_clock_1(t), t))*l_clock_1*m_clock_1*cos(theta_clock_1(t))/(M_clock_1 + m_clock_1))*l_clock_1*m_clock_1*cos(theta_clock_1(t))/((l_clock_1^2*m_clock_1^2*cos(theta_clock_1(t))^2/(M_clock_1 + m_clock_1) - l_clock_1^2*m_clock_1)*(M_clock_1 + m_clock_1)) + (l_clock_1*m_clock_1*sin(theta_clock_1(t))*diff(theta_clock_1(t), t)^2 + k_spring_0*(x_clock_0(t) - x_clock_1(t)) - 2*rho_clock_1*diff(x_clock_1(t), t))/(M_clock_1 + m_clock_1)