Ballistic pendulum
A bullet that sticks in the block of a pendulum at rest, which then swings up.
- Name:
ballistic_pendulum - Field:
physics/mechanics - Formalism: ordinary differential equations
- References: Ryzhov et al. 2021
Factory
@register_model(
field='physics/mechanics',
tags=('mechanics', 'events', 'impacts', 'momentum', 'hybrid'),
references=('Ryzhov et al. 2021',),
defaults={'until': 3.0, 'dt': 0.01},
)
def ballistic_pendulum(
speed: float = 200.0,
distance: float = 1.0,
*,
m: float = 0.01,
M: float = 1.0,
l: float = 1.0,
g: float = 9.81,
) -> System:
"""Return a bullet that sticks in the block of a pendulum at rest, which then swings up.
The bullet flies horizontally towards the block, and the impact is an event. Momentum is conserved in the
perfectly inelastic collision, so the pendulum starts with the speed `m v / (m + M)` and rises to the height
`(m v / (m + M))² / (2 g)`. The bullet position stops at the block after the impact.
Args:
speed:
The speed of the bullet.
distance:
The initial distance of the bullet from the block.
m:
The mass of the bullet.
M:
The mass of the block.
l:
The length of the rod.
g:
The gravitational acceleration.
Returns:
system:
The model, with the bullet position `pendulum.x` and the angle `pendulum.theta`, both of order 2.
Examples:
>>> from skcomplex.models.physics.mechanics import ballistic_pendulum
>>> ballistic_pendulum(distance=2.0).elements[0].x
-2.0
"""
swing = EulerLagrangeInteraction(
'swing',
target='pendulum',
states='theta',
lagrangian=lambda s: s.target.l * derivative(s.target.theta, 1) ** 2 / 2 + s.target.g * cos(s.target.theta),
)
impact = CrossingInteraction(
'impact',
target='pendulum',
guard=lambda s: s.target.x,
reset={
'x_d1': lambda s: 0,
'theta_d1': lambda s: s.target.m * s.target.x_d1 / ((s.target.m + s.target.M) * s.target.l),
},
direction=1,
)
return System(
elements=[BallisticPendulumElement('pendulum', x=-distance, x_d1=speed, m=m, M=M, l=l, g=g)],
interactions=[swing, impact],
)
Components
class BallisticPendulumElement(Element):
"""Bullet of mass `m` at the position `x` flying towards the block of mass `M` of a pendulum with the angle `theta`.
The position and the angle have order 2. The bullet flies freely, the pendulum swings, and the impact transfers
the momentum of the bullet to the pendulum.
Args:
label:
The label of the bullet and pendulum.
visual:
How the bullet and pendulum are drawn, or `None` for none.
x:
The initial position of the bullet, negative before the block.
x_d1:
The initial velocity of the bullet.
theta:
The initial angle of the pendulum.
theta_d1:
The initial angular velocity of the pendulum.
m:
The mass of the bullet.
M:
The mass of the block.
l:
The length of the rod.
g:
The gravitational acceleration.
Examples:
>>> from skcomplex.models.physics.mechanics import BallisticPendulumElement
>>> BallisticPendulumElement('pendulum', M=2.0).M
2.0
"""
def __init__(
self,
label: str,
visual: VisualSpec | None = None,
*,
x: Annotated[float, StateSpec(order=2)] = -1.0,
x_d1: float = 0.0,
theta: Annotated[float, StateSpec(period=2 * pi, order=2)] = 0.0,
theta_d1: float = 0.0,
m: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 0.01,
M: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 1.0,
l: 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.x = x
self.x_d1 = x_d1
self.theta = theta
self.theta_d1 = theta_d1
self.m = m
self.M = M
self.l = l
self.g = g
class EulerLagrangeInteraction(Interaction):
"""Equations of motion `d/dt dL/dq' - dL/dq = Q - dR/dq'` of the coordinates `q` of a Lagrangian `L(ns)`.
The coordinates are states of order 2, and the Lagrangian, the Rayleigh dissipation function `R`, and the
generalized forcing `Q` read them and their first derivatives `derivative(q, 1)`. The interaction returns one
equation law per coordinate, which the compiler solves jointly for the second derivatives, so a Lagrangian with a
mass matrix, such as a double pendulum's, needs no inversion by hand. The equations of one coordinate from several
interactions add, so their Lagrangians, dissipations, and forcings sum: a coupling, such as a spring between two
pendulums, is an interaction of its own whose Lagrangian holds only the coupling energy. The coordinates are the
states of every role that binds one element or group, `target` by default, and every component of a vector state
is a coordinate.
A Lagrangian, dissipation, or forcing function takes the namespace of the interaction, with the namespaces of the
roles, the time `t`, and the symbols of the parameters a subclass declares. The forcing returns one generalized
force per coordinate, in the order of the roles and of the states, or one expression for a single coordinate.
Args:
label:
The label of the interaction, unique among the parts of its system.
target:
The path of the element or group whose states are the coordinates.
states:
The name of the state of order 2 that is a coordinate, or the names of several states.
lagrangian:
A function of the namespace that returns the Lagrangian.
dissipation:
A function of the namespace that returns the Rayleigh dissipation function, or `None` for none.
forcing:
A function of the namespace that returns the generalized forces, or `None` for none.
Raises:
ValidationError:
A state is not a name, or a function is not callable.
Examples:
>>> from typing import Annotated
>>> from sage.all__sagemath_symbolics import cos
>>> from skcomplex.core.elements import Element
>>> from skcomplex.core.interactions import EulerLagrangeInteraction
>>> from skcomplex.core.laws.utils import derivative
>>> from skcomplex.core.system import System
>>> from skcomplex.engine.runtime import Simulation
>>> from skcomplex.specs import StateSpec
>>> class AngleElement(Element):
... def __init__(self, label, q: Annotated[float, StateSpec(order=2)] = 0.1, q_d1: float = 0.0):
... super().__init__(label)
... self.q = q
... self.q_d1 = q_d1
>>> swing = EulerLagrangeInteraction(
... 'swing',
... target='bob',
... states='q',
... lagrangian=lambda s: derivative(s.target.q, 1) ** 2 / 2 + cos(s.target.q),
... )
>>> system = System(elements=[AngleElement('bob')], interactions=[swing])
>>> trajectory = Simulation(system, dt=0.01).check().prepare().run(until=3.2).trajectory_
>>> round(float(trajectory.values('bob.q')[-1]), 2)
-0.1
"""
def __init__(
self,
label: str,
*,
target: Annotated[str, RoleSpec()],
states: str | Sequence[str] = 'x',
lagrangian: Callable[[SimpleNamespace], object] | None = None,
dissipation: Callable[[SimpleNamespace], object] | None = None,
forcing: Callable[[SimpleNamespace], object] | None = None,
) -> None:
super().__init__(label)
self.target = target
self.states = states
self.lagrangian = lagrangian
self.dissipation = dissipation
self.forcing = forcing
def _get_names(self) -> tuple[str, ...]:
"""Return the names of the states that are coordinates."""
return (self.states,) if isinstance(self.states, str) else tuple(self.states)
def _check(self) -> None:
"""Raise an error for a state that is not a name, or a function that is not callable.
Raises:
ValidationError:
The type declares states, a state is not a name, or a function is not callable.
"""
super()._check()
names = self._get_names()
if not names or not all(isinstance(name, str) and name.isidentifier() for name in names):
error_msg = f'{type(self).__name__} {self.label!r} reads the states {self.states!r}. Give names of states.'
raise ValidationError(error_msg)
for name in ('lagrangian', 'dissipation', 'forcing'):
value = getattr(self, name)
if value is not None and not callable(value):
error_msg = f'The {name} of {self.label!r} is {value!r}. Give a function of the namespace.'
raise ValidationError(error_msg)
def _get_lagrangian(self, namespace: SimpleNamespace) -> object:
"""Return the Lagrangian, the hook a subclass overrides to derive it from its parameters."""
if self.lagrangian is None:
error_msg = f'{type(self).__name__} {self.label!r} has no Lagrangian. Give a function of the namespace.'
raise ValidationError(error_msg)
return self.lagrangian(namespace)
def _get_coordinates(self, namespace: SimpleNamespace) -> list[Expression]:
"""Return the coordinates, every component of the states of every role that binds one path."""
coordinates = []
for name, spec in self._get_specs().items():
if not isinstance(spec, RoleSpec) or spec.variadic or spec.noise:
continue
for state in self._get_names():
value = getattr(getattr(namespace, name), state)
coordinates += list(value) if isinstance(value, FreeModuleElement) else [value]
return coordinates
def _get_forces(self, namespace: SimpleNamespace, count: int) -> list[Expression]:
"""Return the generalized force of every coordinate, zero without a forcing."""
if self.forcing is None:
return [SR(0)] * count
forces = self.forcing(namespace)
forces = list(forces) if isinstance(forces, FreeModuleElement | list | tuple) else [forces]
if len(forces) != count:
error_msg = f'The forcing of {self.label!r} returns {len(forces)} forces for {count} coordinates.'
raise ValidationError(error_msg)
return [SR(force) for force in forces]
def _get_laws(self, namespace: SimpleNamespace) -> list[BasePart]:
"""Return the Euler-Lagrange equation of every coordinate.
Args:
namespace:
The namespace of the symbols the laws can use, with the namespaces of the roles.
Returns:
laws:
The equation laws of the coordinates.
"""
coordinates = self._get_coordinates(namespace)
velocities = {derivative(q, 1): SR.var(f'_velocity_{index}') for index, q in enumerate(coordinates)}
back = {symbol: item for item, symbol in velocities.items()}
slots = list(velocities.values())
lagrangian = SR(self._get_lagrangian(namespace)).subs(velocities)
dissipation = SR(0 if self.dissipation is None else self.dissipation(namespace)).subs(velocities)
forces = self._get_forces(namespace, len(coordinates))
time = SR(namespace.t)
laws: list[BasePart] = []
for q, slot, force in zip(coordinates, slots, forces, strict=True):
momentum = lagrangian.diff(slot)
change = momentum.diff(time) + sum(
momentum.diff(other) * other_slot + momentum.diff(other_slot) * derivative(other, 2)
for other, other_slot in zip(coordinates, slots, strict=True)
)
lhs = (change - lagrangian.diff(q)).subs(back)
rhs = (force - dissipation.diff(slot)).subs(back)
laws.append(EquationLaw(lhs == rhs, target=q))
return laws
class CrossingInteraction(Interaction):
"""Event at the crossing of zero by a guard, which sets states of the target to new values or ends the run.
The guard and the new values are functions of the namespace, which holds `namespace.target`, the states and
parameters of the target, and the time `namespace.t`. A guard that reads the members of a group gives every member
its own event, which sets the states of that member alone. Solvers locate the crossing within their tolerance, and
the trajectory records every crossing among its events, under the label of the interaction and with the member.
Args:
label:
The label of the interaction, unique among the parts of its system.
target:
The path of the element or group whose states the guard reads and the event sets.
guard:
The function `guard(namespace)` of the expression whose crossing of zero triggers the event.
reset:
A mapping from the name of a state of the target to the function `f(namespace)` of its new value, applied
all at once, or `None` to keep the states.
direction:
`1` for crossings from negative to positive, `-1` for the opposite, and `0` for both.
terminal:
Whether the crossing ends the run.
Raises:
ValidationError:
The direction is not -1, 0, or 1, or the guard or a new value is not a function.
Examples:
>>> from typing import Annotated
>>> from skcomplex.core.elements import Element
>>> from skcomplex.core.interactions import CrossingInteraction
>>> from skcomplex.core.laws import RateLaw
>>> from skcomplex.core.system import System
>>> from skcomplex.engine.runtime import Simulation
>>> from skcomplex.specs import StateSpec
>>> class ClockElement(Element):
... def __init__(self, label, x: Annotated[float, StateSpec()] = 0.0):
... super().__init__(label)
... self.x = x
...
... def _get_laws(self, s):
... return [RateLaw(s.x, 1.0)]
>>> wrap = CrossingInteraction(
... 'wrap', target='clock', guard=lambda s: s.target.x - 1, reset={'x': lambda s: 0}, direction=1
... )
>>> system = System(elements=[ClockElement('clock')], interactions=[wrap])
>>> trajectory = Simulation(system, dt=0.5).check().prepare().run(until=2.2).trajectory_
>>> [round(event.time, 6) for event in trajectory.events]
[1.0, 2.0]
"""
def __init__(
self,
label: str,
*,
target: Annotated[str, RoleSpec()],
guard: Callable[[SimpleNamespace], object],
reset: Mapping[str, Callable[[SimpleNamespace], object]] | None = None,
direction: int = 0,
terminal: bool = False,
) -> None:
super().__init__(label)
self.target = target
self.guard = guard
self.reset = reset
self.direction = direction
self.terminal = terminal
def _check(self) -> None:
"""Raise an error for a direction that is not -1, 0, or 1, or a guard or a new value that is not a function.
Raises:
ValidationError:
The direction is not -1, 0, or 1, or the guard or a new value is not a function.
"""
super()._check()
if self.direction not in DIRECTIONS:
error_msg = f'The direction of crossing {self.label!r} is {self.direction!r}. Use -1, 0, or 1.'
raise ValidationError(error_msg)
if not callable(self.guard) or not all(callable(function) for function in (self.reset or {}).values()):
error_msg = f'The guard and the new values of crossing {self.label!r} must be functions of the namespace.'
raise ValidationError(error_msg)
@staticmethod
def _get_law_value(value: object) -> object:
"""Return what a function of the namespace returned: a function as it is, or a symbolic expression."""
return value if BaseLaw._is_function(value) else SR(value)
def _get_laws(self, namespace: SimpleNamespace) -> list[BaseLaw]:
"""Return the crossing law of the guard with the new values of the states.
Args:
namespace:
The namespace of the symbols the laws can use, with `namespace.target`.
Returns:
laws:
The crossing law.
"""
reset = {
getattr(self._get_member(namespace, 'target'), name): self._get_law_value(f(namespace))
for name, f in (self.reset or {}).items()
}
return [
CrossingLaw(
self._get_law_value(self.guard(namespace)),
reset or None,
self.direction,
self.terminal,
label=self.label,
),
]
Equations
pendulum.x: diff(x_pendulum(t), t, t) == 0
pendulum.theta: diff(theta_pendulum(t), t, t) == -g_pendulum*sin(theta_pendulum(t))/l_pendulum