Dynamical billiards of Sinai, mushroom, and open three-disk 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; Gaspard and Rice 1989; Datseris and Parlitz 2022
Select this system
from skcomplex.systems import create_simulation, create_system, describe
system = create_system('billiards')
entry = describe('billiards')
trajectory = create_simulation('billiards').prepare().run().trajectory_
Factory
@register_system(
field='physics/nonlinear_oscillators',
tags=('events', 'billiard', 'chaos', 'mixed phase space'),
references=('Sinai 1970', 'Bunimovich 2001', 'Gaspard and Rice 1989', '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,
spacing: float = 2.4,
) -> 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. The three-disk table is open: three discs of radius `disc_radius` sit at the corners of an equilateral
triangle of side `spacing` centered at the origin, with one disc below the origin. A particle that comes from far
away is scattered chaotically between the discs and leaves again, and the number of its collisions and its outgoing
angle depend on its incoming direction through a fractal scattering function. Every bounce is an event, and the
speed is conserved.
Args:
shape:
`'sinai'`, `'mushroom'`, or `'three_disk'`.
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 and of the discs of the three-disk 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.
spacing:
The distance between the centers of the discs of the three-disk table.
Returns:
system:
The system, with the position `particle.x`, `particle.y` and the velocity `particle.vx`, `particle.vy`.
Raises:
ValidationError:
The shape is unknown.
Examples:
>>> from skcomplex.systems.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)
if shape == 'three_disk':
table = ThreeDiskBilliardElement('particle', x=x, y=y, vx=vx, vy=vy, spacing=spacing, radius=disc_radius)
walls = [_disc(f'disc_{index}', angle) for index, angle in enumerate(ANGLES)]
return System(elements=[table], interactions=walls)
error_msg = f'Unknown billiard shape {shape!r}. Use one of {", ".join(SHAPES)}.'
raise ValidationError(error_msg)
Components
Elements:
SinaiBilliardElement, fromElement
Interactions:
ReflectionInteraction, fromCrossingInteraction
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.systems.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.components.elements import PointElement
>>> from skcomplex.core.components.spaces import ContinuousSpace
>>> from skcomplex.core.components.system import System
>>> from skcomplex.core.engine.runtime import Simulation
>>> from skcomplex.systems.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).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