Bouncing ball
A ball that falls under gravity and bounces off the floor, losing speed by a coefficient of restitution.
- Name:
bouncing_ball - Field:
physics/mechanics - Formalism: ordinary differential equations
- References: Strogatz 2015; Hairer, Nørsett, and Wanner 1993
Factory
@register_model(
field='physics/mechanics',
tags=('mechanics', 'events', 'impacts'),
references=('Strogatz 2015', 'Hairer, Nørsett, and Wanner 1993'),
defaults={'until': 3.0, 'dt': 0.01},
)
def bouncing_ball(h: float = 1.0, v: float = 0.0, *, g: float = 9.81, e: float = 0.8) -> System:
"""Return a ball that falls under gravity and bounces off the floor, losing speed by a coefficient of restitution.
Every impact is an event that reverses the velocity and scales it by `e`. The impacts accumulate at the Zeno time
`sqrt(2 h / g) (1 + e) / (1 - e)` for a ball dropped from rest, so runs should end before it.
Args:
h:
The initial height.
v:
The initial vertical velocity.
g:
The gravitational acceleration.
e:
The coefficient of restitution.
Returns:
system:
The model, with the height `ball.h` of order 2 and its velocity `ball.h_d1`.
Examples:
>>> from skcomplex.models.physics.mechanics import bouncing_ball
>>> from skcomplex.engine.runtime import Simulation
>>> trajectory = Simulation(bouncing_ball(), dt=0.01).check().prepare().run(until=1.0).trajectory_
>>> round(trajectory.events[0].time, 6)
0.451524
"""
return System(
elements=[BouncingBallElement('ball', h=h, h_d1=v, g=g)],
interactions=[
RateInteraction('gravity', target='ball', state='h', rate=lambda s: -s.target.g),
ReflectionInteraction(
'impact',
target='ball',
guard=lambda s: s.target.h,
positions=('h',),
velocities=('h_d1',),
direction=-1,
restitution=e,
),
],
)
Components
class BouncingBallElement(Element):
"""Ball at the height `h` above a floor, falling under the gravitational acceleration `g`.
The height has order 2 and its first derivative `h_d1` is the vertical velocity.
Args:
label:
The label of the ball.
visual:
How the ball is drawn, or `None` for none.
h:
The initial height.
h_d1:
The initial vertical velocity.
g:
The gravitational acceleration.
Examples:
>>> from skcomplex.models.physics.mechanics import BouncingBallElement
>>> BouncingBallElement('ball', h=2.0).h
2.0
"""
def __init__(
self,
label: str,
visual: VisualSpec | None = None,
*,
h: Annotated[float, StateSpec(order=2)] = 1.0,
h_d1: float = 0.0,
g: Annotated[float, ParameterSpec()] = 9.81,
) -> None:
super().__init__(label, visual=visual)
self.h = h
self.h_d1 = h_d1
self.g = g
class RateInteraction(Interaction):
"""Rate `rate(ns)` and noise `noise(ns) dW` added to the highest time derivative of a state of a target.
On a state of order 1 the rate is its first derivative, and on a state of order `k` it adds to the `k`-th
derivative, so rates on one state superpose. With a divisor `d` the rate and the noise are divided by `d`, the name
of a state or parameter of the target, such as an inertia, a friction, or a capacitance, which backends can read.
A rate, a noise, or a divisor function takes the namespace of the interaction, with the namespaces `target` and
`sources` of the bound elements, the time `t`, and the symbols of the parameters a subclass declares. A number
stands for a constant, and a number on a vector state applies to every component. The noise is a Wiener process per
component and member by default, a process of a `WienerProcessSpec` to share or correlate it across states and
members, or a fractional Brownian motion of a Hurst index.
Args:
label:
The label of the interaction, unique among the parts of its system.
target:
The path of the element or group whose state the rate changes.
state:
The name of the state of the target that the rate changes.
rate:
A function of the namespace that returns the rate, or a number.
sources:
The paths of the elements or groups the rate and the noise read but do not change.
noise:
A function of the namespace that returns the amplitude of the noise, a number, or `None` for no noise.
divisor:
The name of the state or parameter of the target that divides the rate and the noise, a function of the
namespace of the target that returns the divisor, such as a moment of inertia `m l^2`, or `None`.
process:
The Wiener processes that drive the noise, or `None` for one independent process per component and member.
index:
The process of `process` that drives the noise, or `None` for every process, one per component of a vector
state.
convention:
The stochastic calculus of the noise, `'ito'`, `'stratonovich'`, or `'anti-ito'`.
hurst:
The Hurst index in `(0, 1)` of a fractional Brownian motion that drives the noise, a number or a function
of the namespace that returns a parameter, or `None` for a Wiener process.
conserved:
Whether the noise is the divergence of a random flux on the links of the lattice of a field, so it keeps the
total of the state.
Raises:
ValidationError:
The state is not a name, a function is neither callable nor a number, the convention is unknown, the Hurst
index is outside `(0, 1)`, or the noise options cannot combine.
Examples:
>>> from skcomplex.core.elements import CompartmentElement
>>> from skcomplex.core.interactions import RateInteraction
>>> from skcomplex.core.system import System
>>> from skcomplex.engine.runtime import Simulation
>>> decay = RateInteraction('decay', target='stock', state='n', rate=lambda s: -0.5 * s.target.n)
>>> system = System(elements=[CompartmentElement('stock', n=8.0)], interactions=[decay])
>>> trajectory = Simulation(system, dt=0.01).check().prepare().run(until=2.0).trajectory_
>>> round(float(trajectory.values('stock.n')[-1]), 3)
2.943
"""
def __init__(
self,
label: str,
*,
target: Annotated[str, RoleSpec()],
state: str = 'x',
rate: Callable[[SimpleNamespace], object] | float = 0.0,
sources: Annotated[Sequence[str], RoleSpec(variadic=True)] = (),
noise: Callable[[SimpleNamespace], object] | float | None = None,
divisor: str | Callable[[SimpleNamespace], object] | None = None,
process: Annotated[WienerProcessSpec | None, RoleSpec(noise=True)] = None,
index: int | None = None,
convention: str = 'ito',
hurst: Callable[[SimpleNamespace], object] | float | None = None,
conserved: bool = False,
) -> None:
super().__init__(label)
self.target = target
self.state = state
self.rate = rate
self.sources = sources
self.noise = noise
self.divisor = divisor
self.process = process
self.index = index
self.convention = convention
self.hurst = hurst
self.conserved = conserved
def _check_functions(self) -> None:
"""Raise an error for a state that is not a name, or a function that is neither callable nor a number."""
if not isinstance(self.state, str) or not self.state.isidentifier():
error_msg = (
f'{type(self).__name__} {self.label!r} changes the state {self.state!r}, which is not a name of a '
'state.'
)
raise ValidationError(error_msg)
for name in ('rate', 'noise'):
value = getattr(self, name)
if not (value is None or callable(value) or isinstance(value, int | float)):
error_msg = f'The {name} of {self.label!r} is {value!r}. Give a function of the namespace or a number.'
raise ValidationError(error_msg)
if not (
self.divisor is None
or callable(self.divisor)
or (isinstance(self.divisor, str) and self.divisor.isidentifier())
):
error_msg = (
f'The divisor of {self.label!r} is {self.divisor!r}. Give the name of a state or parameter, or a '
'function of the namespace of the target.'
)
raise ValidationError(error_msg)
def _check_noise(self) -> None:
"""Raise an error for an unknown convention, a Hurst index outside `(0, 1)`, or noise options that conflict."""
if self.convention not in CONVENTIONS:
error_msg = f'Unknown noise convention {self.convention!r}. Use one of {", ".join(CONVENTIONS)}.'
raise ValidationError(error_msg)
if not (
self.hurst is None or callable(self.hurst) or (isinstance(self.hurst, int | float) and 0 < self.hurst < 1)
):
error_msg = f'The Hurst index of {self.label!r} is {self.hurst!r}. Give a number in (0, 1) or a function.'
raise ValidationError(error_msg)
if self.index is not None and self.process is None:
error_msg = f'{self.label!r} picks the process {self.index!r} but binds no processes. Give `process`.'
raise ValidationError(error_msg)
options = [
name
for name, given in (('hurst', self.hurst), ('process', self.process), ('conserved', self.conserved))
if given
]
if len(options) > 1:
error_msg = f'The noise of {self.label!r} takes {" and ".join(options)}. Give at most one of them.'
raise ValidationError(error_msg)
def _check(self) -> None:
"""Raise an error for a state, functions, or noise options the interaction cannot take.
Raises:
ValidationError:
The type declares states, the state is not a name, a function is neither callable nor a number, the
convention is unknown, the Hurst index is outside `(0, 1)`, or the noise options cannot combine.
"""
super()._check()
self._check_functions()
self._check_noise()
def _evaluate(self, function: Callable[[SimpleNamespace], object] | float, namespace: SimpleNamespace) -> object:
"""Return the value of a function of the namespace, or a number as it is."""
return function(namespace) if callable(function) else function
def _broadcast(self, value: object, state: object) -> object:
"""Return a value with one component per component of a vector state, repeating a number."""
if isinstance(state, FreeModuleElement) and not isinstance(value, FreeModuleElement | list | tuple):
return vector(SR, [value] * len(state))
if isinstance(value, list | tuple):
return vector(SR, value)
return value
def _divide(self, value: object, namespace: SimpleNamespace) -> object:
"""Return a value divided by the divisor of the target, or the value without a divisor."""
if self.divisor is None:
return value
divisor = (
self.divisor(self._get_member(namespace, 'target'))
if callable(self.divisor)
else getattr(self._get_member(namespace, 'target'), self.divisor)
)
return operator.truediv(value, divisor)
def _get_rate(self, namespace: SimpleNamespace) -> object:
"""Return the rate before the divisor, the hook a subclass overrides to derive it from its parameters."""
return self._evaluate(self.rate, namespace)
def _get_noise(self, namespace: SimpleNamespace) -> object:
"""Return the noise amplitude before the divisor, the hook a subclass overrides, or `None` for no noise."""
return None if self.noise is None else self._evaluate(self.noise, namespace)
def _get_noise_law(self, state: object, amplitude: object, namespace: SimpleNamespace) -> BasePart:
"""Return the noise law of an amplitude, by the process, the convention, and the Hurst index."""
if self.hurst is not None:
hurst = self._evaluate(self.hurst, namespace)
return FractionalNoiseLaw(state, amplitude, hurst=hurst, convention=self.convention)
process = None
if self.process is not None:
handle = self._get_member(namespace, 'process').dW
process = handle if self.index is None else handle[self.index]
return NoiseLaw(state, amplitude, convention=self.convention, process=process, conserved=self.conserved)
def _get_laws(self, namespace: SimpleNamespace) -> list[BasePart]:
"""Return the rate law of the state and its noise law, both divided by the divisor.
Args:
namespace:
The namespace of the symbols the laws can use, with the namespaces of the target and the sources.
Returns:
laws:
The rate law, and the noise law when the interaction has noise.
"""
state = getattr(self._get_member(namespace, 'target'), self.state)
rate = self._divide(self._broadcast(self._get_rate(namespace), state), namespace)
laws: list[BasePart] = [RateLaw(state, rate)]
noise = self._get_noise(namespace)
if noise is not None:
laws.append(self._get_noise_law(state, self._divide(noise, namespace), namespace))
return laws
class ReflectionInteraction(CrossingInteraction):
"""Reflection of a moving point off a wall `g(x) = 0`, with a coefficient of restitution `e`.
When the guard `g` crosses zero, the velocity `v` of the point becomes `v - (1 + e) (v . n) n / |n|**2`, where `n`
is the gradient of `g` with respect to the position. The normal velocity is reversed and scaled by `e`, and the
tangential velocity is kept, so `e = 1` is a specular, elastic reflection. The position and the velocity are given
by the names of their states, scalar or vector, whose components pair up in order. A `region` restricts the wall to
the part of the surface where it holds, and the point crosses the rest of the surface freely.
Args:
label:
The label of the interaction.
target:
The path of the element that reflects.
guard:
The function `guard(namespace)` of the wall, an expression of the position that is zero on the wall.
positions:
The names of the states that hold the position.
velocities:
The names of the states that hold the velocity, the first derivatives of the position.
direction:
`1` for crossings from negative to positive, `-1` for the opposite, and `0` for both.
restitution:
The coefficient of restitution, the ratio of the normal speeds after and before the reflection.
region:
The function `region(namespace)` of the condition that holds on the part of the surface that is a wall, or
`None` for a wall on the whole surface.
Raises:
ValidationError:
The positions and the velocities have different numbers of names, or a function is not callable.
Examples:
>>> from skcomplex.core.elements import PointElement
>>> from skcomplex.core.spaces import ContinuousSpace
>>> from skcomplex.core.system import System
>>> from skcomplex.engine.runtime import Simulation
>>> from skcomplex.models.physics.common.mechanics import ReflectionInteraction
>>> wall = ReflectionInteraction('wall', target='ball', guard=lambda s: s.target.x[0] - 1, direction=1)
>>> ball = PointElement('ball', x=[0.0, 0.0], x_d1=[1.0, 1.0])
>>> system = System(elements=[ball], interactions=[wall], space=ContinuousSpace(n_dim=2))
>>> trajectory = Simulation(system, dt=0.5).check().prepare().run(until=2.0).trajectory_
>>> trajectory.values('ball.x')[-1].tolist()
[0.0, 2.0]
"""
def __init__(
self,
label: str,
*,
target: Annotated[str, RoleSpec()],
guard: Callable[[SimpleNamespace], object],
positions: Sequence[str] = ('x',),
velocities: Sequence[str] = ('x_d1',),
direction: int = 0,
restitution: Annotated[float, ParameterSpec(domain=RealSet.closed(0, 1))] = 1.0,
region: Callable[[SimpleNamespace], object] | None = None,
) -> None:
super().__init__(label, target=target, guard=guard, direction=direction)
self.positions = positions
self.velocities = velocities
self.restitution = restitution
self.region = region
def _get_components(self, member: SimpleNamespace, names: Sequence[str]) -> list[Expression]:
"""Return the components of the named states of a member, every component of a vector state."""
components: list[Expression] = []
for name in names:
value = getattr(member, name)
components += list(value) if isinstance(value, FreeModuleElement) else [SR(value)]
return components
def _check(self) -> None:
"""Raise an error for positions and velocities of different numbers of names, or a function not callable.
Raises:
ValidationError:
The positions and the velocities have different numbers of names, or a function is not callable.
"""
super()._check()
if self.region is not None and not callable(self.region):
error_msg = f'The region of reflection {self.label!r} must be a function of the namespace.'
raise ValidationError(error_msg)
if len(self.positions) != len(self.velocities):
error_msg = (
f'Reflection {self.label!r} pairs the positions {self.positions!r} with the velocities '
f'{self.velocities!r}. Give one velocity per position.'
)
raise ValidationError(error_msg)
def _get_laws(self, namespace: SimpleNamespace) -> list[BaseLaw]:
"""Return the crossing law of the wall, which reflects the velocity about the normal of the wall.
Args:
namespace:
The namespace of the symbols the laws can use, with `namespace.target`.
Returns:
laws:
The crossing law.
"""
guard = SR(self.guard(namespace))
positions = self._get_components(namespace.target, self.positions)
velocities = self._get_components(namespace.target, self.velocities)
normal = [guard.diff(position) for position in positions]
speed = sum((v * n for v, n in zip(velocities, normal, strict=True)), SR(0))
scale = (1 + namespace.restitution) * speed / sum((n**2 for n in normal), SR(0))
reset = {v: v - scale * n for v, n in zip(velocities, normal, strict=True)}
if self.region is not None:
region = self.region(namespace)
reset = {v: where(region, value, v) for v, value in reset.items()}
return [CrossingLaw(guard, reset, self.direction, self.terminal, label=self.label)]
Equations
ball.h: diff(h_ball(t), t, t) == -g_ball