Coupled pendulums
Two pendulums whose bobs are joined by a spring, which pass their energy back and forth.
- Name:
coupled_pendulums - Field:
physics/mechanics - Formalism: ordinary differential equations
- References: Ryzhov et al. 2021
Factory
@register_model(
field='physics/mechanics',
tags=('mechanics', 'ode', 'second order', 'oscillation', 'normal modes'),
references=('Ryzhov et al. 2021',),
defaults={'until': 40.0, 'dt': 0.05},
)
def coupled_pendulums(
theta1: float = 0.2,
theta2: float = 0.0,
*,
k: float = 0.5,
m1: float = 1.0,
m2: float = 1.0,
l1: float = 1.0,
l2: float = 1.0,
g: float = 9.81,
) -> System:
"""Return two pendulums whose bobs are joined by a spring, which pass their energy back and forth.
The spring stretches with the horizontal distance `l2 sin θ2 - l1 sin θ1` of the bobs and adds the torque
`k l1 (l2 sin θ2 - l1 sin θ1) cos θ1` to the first pendulum and the opposite one to the second. Identical pendulums
have the normal modes `sqrt(g / l)` in phase and `sqrt(g / l + 2 k / m)` in antiphase.
Args:
theta1:
The initial angle of the first pendulum.
theta2:
The initial angle of the second pendulum.
k:
The stiffness of the spring.
m1:
The mass of the first bob.
m2:
The mass of the second bob.
l1:
The length of the first pendulum.
l2:
The length of the second pendulum.
g:
The gravitational acceleration.
Returns:
system:
The model, with the pendulums `left` and `right` and their angles `left.theta` and `right.theta`.
Examples:
>>> from skcomplex.models.physics.mechanics import coupled_pendulums
>>> coupled_pendulums(k=2.0).interactions[0].k
2.0
"""
return System(
elements=[
PendulumElement('left', theta=theta1, m=m1, l=l1, g=g),
PendulumElement('right', theta=theta2, m=m2, l=l2, g=g),
],
interactions=[
PendulumSpringInteraction('spring', first='left', second='right', k=k),
PendulumInteraction('left_swing', target='left'),
PendulumInteraction('right_swing', target='right'),
],
)
Components
class PendulumElement(Element):
"""Pendulum, a bob of mass `m` on a rigid rod of length `l` that swings by the angle `theta` under gravity `g`.
The angle has order 2 and its first derivative `theta_d1` is the angular velocity. The motion comes from an
interaction such as `PendulumInteraction`.
Args:
label:
The label of the pendulum.
visual:
How the pendulum is drawn, or `None` for its class visual.
theta:
The initial angle from the vertical, in radians.
theta_d1:
The initial angular velocity.
m:
The mass of the bob.
g:
The gravitational acceleration.
l:
The length of the rod.
Examples:
>>> from skcomplex.models.physics.mechanics import PendulumElement
>>> PendulumElement('bob', theta=0.2, l=2.0).l
2.0
"""
visual = VisualSpec('bodies', {'positions': _pendulum_bob, 'pivot': (0.0, 0.0), 'rods': True})
def __init__(
self,
label: str,
visual: VisualSpec | None = None,
*,
theta: Annotated[float, StateSpec(period=2 * pi, order=2)] = 0.5,
theta_d1: float = 0.0,
m: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 1.0,
g: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 9.81,
l: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 1.0,
) -> None:
super().__init__(label, visual=visual)
self.theta = theta
self.theta_d1 = theta_d1
self.m = m
self.g = g
self.l = l
class PendulumSpringInteraction(EulerLagrangeInteraction):
"""Spring between the bobs of two pendulums, which stretches with the horizontal distance of the bobs.
The spring of stiffness `k` adds the Lagrangian `-k (l2 sin θ2 - l1 sin θ1)² / 2` to the Euler-Lagrange equations
of the two angles, so its torque on every pendulum divides by the moment of inertia `m l²` of the pendulum's own
Lagrangian.
Args:
label:
The label of the interaction.
first:
The path of the first pendulum, with the angle `theta`, the mass `m`, and the length `l`.
second:
The path of the second pendulum.
k:
The stiffness of the spring.
Examples:
>>> from skcomplex.models.physics.mechanics import PendulumSpringInteraction
>>> PendulumSpringInteraction('spring', first='left', second='right', 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='theta')
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.
"""
first, second = namespace.first, namespace.second
stretch = second.l * sin(second.theta) - first.l * sin(first.theta)
return -namespace.k * stretch**2 / 2
class PendulumInteraction(EulerLagrangeInteraction):
"""Motion of a pendulum, the Euler-Lagrange equation of `L = m l² θ'² / 2 + m g l cos θ`.
The equation is Newton's law `m l² θ'' = -m g l sin θ` for the angle `θ` of the target, whose mass `m`, length
`l`, and gravity `g` are parameters of the target.
Args:
label:
The label of the interaction.
target:
The path of the pendulum.
Examples:
>>> from skcomplex.models.physics.mechanics import PendulumInteraction
>>> PendulumInteraction('swing', target='bob').states
'theta'
"""
def __init__(self, label: str, *, target: Annotated[str, RoleSpec()]) -> None:
super().__init__(label, target=target, states='theta')
def _get_lagrangian(self, namespace: SimpleNamespace) -> Expression:
"""Return the Lagrangian `m l² θ'² / 2 + m g l cos θ` of the pendulum.
Args:
namespace:
The namespace of the symbols the Lagrangian can use.
Returns:
lagrangian:
The Lagrangian.
"""
bob = namespace.target
return bob.m * bob.l**2 * derivative(bob.theta, 1) ** 2 / 2 + bob.m * bob.g * bob.l * cos(bob.theta)
Equations
left.theta: diff(theta_left(t), t, t) == -((l_left*sin(theta_left(t)) - l_right*sin(theta_right(t)))*k_spring*l_left*cos(theta_left(t)) + g_left*l_left*m_left*sin(theta_left(t)))/(l_left^2*m_left)
right.theta: diff(theta_right(t), t, t) == ((l_left*sin(theta_left(t)) - l_right*sin(theta_right(t)))*k_spring*l_right*cos(theta_right(t)) - g_right*l_right*m_right*sin(theta_right(t)))/(l_right^2*m_right)