Dynamical billiards of Sinai and mushroom shapes
A particle moving freely on a billiard table and reflecting specularly off its walls.
- Name:
billiards - Field:
physics/nonlinear_oscillators - Formalism: ordinary differential equations
- References: Sinai 1970; Bunimovich 2001; Datseris and Parlitz 2022
Factory
@register_model(
field='physics/nonlinear_oscillators',
tags=('events', 'billiard', 'chaos', 'mixed phase space'),
references=('Sinai 1970', 'Bunimovich 2001', 'Datseris and Parlitz 2022'),
defaults={'until': 20.0, 'dt': 0.05},
)
def billiards(
shape: str = 'sinai',
x: float = 0.1,
y: float = 0.6,
vx: float = 0.6,
vy: float = 0.8,
*,
size: float = 1.0,
disc_radius: float = 0.4,
cap_radius: float = 1.0,
width: float = 0.4,
height: float = 1.0,
) -> System:
"""Return a particle moving freely on a billiard table and reflecting specularly off its walls.
The Sinai table is the square `[-size, size]²` with a disc of radius `disc_radius` at its center, and its motion is
chaotic. The mushroom table is the upper half of a disc of radius `cap_radius` on a stem of half width `width` and
height `height`, and its phase space mixes regular orbits that stay in the cap with chaotic orbits that enter the
stem. Every bounce is an event, and the speed is conserved.
Args:
shape:
`'sinai'` or `'mushroom'`.
x:
The initial horizontal position.
y:
The initial vertical position.
vx:
The initial horizontal velocity.
vy:
The initial vertical velocity.
size:
The half side of the Sinai table.
disc_radius:
The radius of the disc of the Sinai table.
cap_radius:
The radius of the cap of the mushroom.
width:
The half width of the stem of the mushroom.
height:
The height of the stem of the mushroom.
Returns:
system:
The model, with the position `particle.x`, `particle.y` and the velocity `particle.vx`, `particle.vy`.
Raises:
ValidationError:
The shape is unknown.
Examples:
>>> from skcomplex.models.physics.nonlinear_oscillators import billiards
>>> type(billiards('mushroom').elements[0]).__name__
'MushroomBilliardElement'
"""
if shape == 'sinai':
table = SinaiBilliardElement('particle', x=x, y=y, vx=vx, vy=vy, size=size, radius=disc_radius)
walls = [
_wall('right_wall', lambda s: s.target.x - s.target.size, 1),
_wall('left_wall', lambda s: s.target.x + s.target.size, -1),
_wall('top_wall', lambda s: s.target.y - s.target.size, 1),
_wall('bottom_wall', lambda s: s.target.y + s.target.size, -1),
_wall('disc', lambda s: s.target.x**2 + s.target.y**2 - s.target.radius**2, -1),
]
return System(elements=[table], interactions=walls)
if shape == 'mushroom':
table = MushroomBilliardElement(
'particle',
x=x,
y=y,
vx=vx,
vy=vy,
radius=cap_radius,
width=width,
height=height,
)
walls = [
_wall('cap', lambda s: s.target.x**2 + s.target.y**2 - s.target.radius**2, 1, lambda s: s.target.y > 0),
_wall('cap_floor', lambda s: s.target.y, -1, lambda s: abs(s.target.x) > s.target.width),
_wall('right_stem_wall', lambda s: s.target.x - s.target.width, 1, lambda s: s.target.y < 0),
_wall('left_stem_wall', lambda s: s.target.x + s.target.width, -1, lambda s: s.target.y < 0),
_wall('stem_floor', lambda s: s.target.y + s.target.height, -1),
]
return System(elements=[table], interactions=walls)
error_msg = f'Unknown billiard shape {shape!r}. Use one of {", ".join(SHAPES)}.'
raise ValidationError(error_msg)
Components
class SinaiBilliardElement(Element):
"""Free particle in a square table of half side `size` with a central disc of radius `radius`, the Sinai billiard.
The particle flies straight, and its specular reflections off the four walls and the disc make its motion ergodic
and chaotic.
Args:
label:
The label of the particle.
visual:
How the particle is drawn.
x:
The initial first coordinate.
y:
The initial second coordinate.
vx:
The initial first velocity.
vy:
The initial second velocity.
size:
The half side of the square.
radius:
The radius of the disc.
Attributes:
default_visual_:
The visual of the particle when it sets none.
Examples:
>>> from skcomplex.models.physics.nonlinear_oscillators import SinaiBilliardElement
>>> SinaiBilliardElement('ball').radius
0.4
"""
def __init__(
self,
label: str,
visual: VisualSpec | None = VisualSpec(
'bodies',
{'positions': ('x', 'y'), 'trail': None, 'outline': _sinai_walls},
),
*,
x: Annotated[float, StateSpec()] = 0.1,
y: Annotated[float, StateSpec()] = 0.6,
vx: Annotated[float, StateSpec()] = 0.6,
vy: Annotated[float, StateSpec()] = 0.8,
size: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 1.0,
radius: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 0.4,
) -> None:
super().__init__(label, visual=visual)
self.x = x
self.y = y
self.vx = vx
self.vy = vy
self.size = size
self.radius = radius
def _get_laws(self, namespace: SimpleNamespace) -> list[BaseLaw]:
"""Return the free flight `x' = vx`, `y' = vy` of the particle.
Args:
namespace:
The namespace of the symbols the laws can use.
Returns:
laws:
The rates of the position.
"""
s = namespace
return [RateLaw(s.x, s.vx), RateLaw(s.y, s.vy)]
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
particle.x: diff(x_particle(t), t) == vx_particle(t)
particle.y: diff(y_particle(t), t) == vy_particle(t)
particle.vx: diff(vx_particle(t), t) == 0
particle.vy: diff(vy_particle(t), t) == 0