physicskit.chaos#

physicskit.chaos hub · History · Examples · API reference

physicskit.chaos: visual analysis and simulation of chaotic dynamical systems and 2D billiards.

class physicskit.chaos.BakersMap(alpha=0.5)[source]#

Bases: DiscreteMap

The (generalized) baker’s map on the unit square: the textbook model of chaos.

The map cuts the unit square [0, 1) x [0, 1) at x = alpha, stretches each piece horizontally back to unit width (contracting it vertically to match), and stacks the two pieces – literally the “stretch, cut, and stack” mechanism used to motivate deterministic chaos:

\[\begin{split}T(x, y) = \begin{cases} (x / \alpha,\ \alpha y) & 0 \le x < \alpha \\ ((x - \alpha) / (1 - \alpha),\ \alpha + (1 - \alpha) y) & \alpha \le x < 1 \end{cases}\end{split}\]

Because each branch is affine with Jacobian determinant exactly 1, the map is area-preserving (unlike the dissipative Henon map) while still being uniformly hyperbolic, ergodic, and mixing: a single long orbit fills the unit square uniformly and densely, and its Lyapunov exponents are known exactly in closed form (see lyapunov_exponents()), making it the standard textbook example for validating numerical chaos estimators. alpha=0.5 (the default) recovers the classic symmetric baker’s map.

Parameters:

alpha (float) – Cut position, in (0, 1); the classic symmetric map has alpha=0.5.

Variables:

alpha (float) – Cut position.

Raises:

ValueError – If alpha does not satisfy 0 < alpha < 1.

Notes

x = 0 is an exact fixed point of the map for every alpha (since 0 / alpha = 0, always in the first branch). Because each branch is expanding, floating-point rounding error is amplified every step, and for many choices of alpha and the initial condition, a long-enough orbit eventually rounds to exactly 0.0 in floating point and gets permanently trapped there – silently producing a degenerate, physically meaningless trajectory rather than raising an error. This typically happens within a few hundred to a few thousand iterations (it depends sensitively on alpha and the initial condition: it may also never happen within a given run). Do not rely on trajectory() orbits of more than a few hundred iterations to be a faithful sample of the invariant measure without first checking the orbit hasn’t collapsed; HenonMap and continuous systems such as Lorenz do not share this failure mode and are better suited to long-trajectory statistics.

dim: int = 2#

State dimension, always 2. State is (x, y).

initial_state()[source]#

Default initial condition (0.5, 0.5).

Return type:

NDArray[double]

Returns:

ndarray of float, shape (2,)

lyapunov_exponents()[source]#

Exact Lyapunov exponents of the map, in closed form.

The baker’s map is piecewise-linear and uniformly hyperbolic: almost every orbit spends a fraction alpha of its time in the expanding-by 1/alpha branch and 1 - alpha in the expanding-by 1/(1 - alpha) branch, so by the ergodic theorem its Lyapunov exponent equals the Shannon entropy of that two-symbol Bernoulli process – no numerical estimation needed. This is exactly what makes the baker’s map useful as a ground truth for validating physicskit.chaos.utils.metrics.lyapunov_exponent_from_divergence() and similar numerical estimators.

Return type:

tuple[float, float]

Returns:

lambda_expanding, lambda_contracting (float) – The two Lyapunov exponents, +h and -h, where h = -alpha*ln(alpha) - (1-alpha)*ln(1-alpha) (= ln(2) for the classic symmetric map, alpha=0.5).

step(state)[source]#

Advance (x, y) by one map iteration.

Parameters:

state (NDArray[double]) – Current state (x, y), each in [0, 1).

Return type:

NDArray[double]

Returns:

ndarray of float, shape (2,) – Next state (x, y).

trajectory(state0=None, n_iter=1000)[source]#

Iterate the map n_iter times starting from state0.

Parameters:
Return type:

NDArray[double]

Returns:

ndarray of float, shape (n_iter + 1, 2) – State at each iteration, including state0 as row 0.

class physicskit.chaos.BilliardSystem[source]#

Bases: ABC

Base class for 2D billiard geometries.

A billiard is described as a closed boundary made of straight segments and circular arcs. Subclasses build this boundary representation once (in _build_boundary) and this base class provides the shared, Numba-accelerated ray-tracing / specular-reflection machinery plus boundary-coordinate (arclength s, sine of the reflection angle sin(phi)) bookkeeping used for Poincare sections.

Variables:

  • _segments (ndarray of float, shape (n_segments, 4)) – Straight boundary walls, each row (x1, y1, x2, y2).

  • _arcs (ndarray of float, shape (n_arcs, 5)) – Circular-arc boundary walls, each row (cx, cy, r, theta1, theta2).

  • _arc_full (ndarray of bool, shape (n_arcs,)) – Whether each arc is a full circle (True) or a bounded arc.

  • _s_offsets_seg (ndarray of float, shape (n_segments,)) – Cumulative boundary arclength at the start of each segment.

  • _s_offsets_arc (ndarray of float, shape (n_arcs,)) – Cumulative boundary arclength at the start of each arc.

  • _perimeter (float) – Total boundary perimeter.

boundary_arrays()[source]#

Return the raw boundary wall arrays.

Return type:

tuple[NDArray[double], NDArray[double], NDArray[bool]]

Returns:

  • segments (ndarray of float, shape (n_segments, 4)) – Straight walls, each row (x1, y1, x2, y2).

  • arcs (ndarray of float, shape (n_arcs, 5)) – Circular-arc walls, each row (cx, cy, r, theta1, theta2).

  • arc_full (ndarray of bool, shape (n_arcs,)) – Whether each arc is a full circle.

abstractmethod boundary_polyline(points_per_arc=200)[source]#

Trace the boundary as one or more closed polylines, suitable for plotting.

Parameters:

points_per_arc (int) – Number of points used to sample each circular arc.

Return type:

NDArray[double]

Returns:

ndarray of float, shape (M, 2) – Points tracing the boundary; disjoint closed components (if any) are separated by a row of NaN. See the concrete implementation (shared by every billiard shape) for details.

perimeter()[source]#

Total boundary perimeter.

Return type:

float

Returns:

float – The billiard boundary’s total arclength.

abstractmethod sample_interior_point()[source]#

Return a point guaranteed to lie in the billiard’s interior.

Return type:

NDArray[double]

Returns:

ndarray of float, shape (2,) – An (x, y) point strictly inside the billiard.

abstractmethod simulate(pos, vel, n_bounces)[source]#

Trace n_bounces specular reflections from an initial ray.

Parameters:

  • pos (ArrayLike) – Initial position (x, y); must lie in the billiard’s interior.

  • vel (ArrayLike) – Initial velocity direction (vx, vy); normalized internally.

  • n_bounces (int) – Number of reflections to simulate.

Return type:

dict[str, NDArray[Any]]

Returns:

dict – Dictionary with keys x, y, vx, vy, s, sin_phi, wall_type, wall_idx, each an array of length n_bounces; see the concrete implementation (shared by every billiard shape) for details.

abstractmethod simulate_many_rays(pos, angles, n_bounces)[source]#

Trace many independent rays from one point, in parallel.

Parameters:

  • pos (ArrayLike) – Shared initial position (x, y) for every ray; must lie in the billiard’s interior.

  • angles (ArrayLike) – Initial launch angle (radians) of each ray.

  • n_bounces (int) – Number of reflections to simulate per ray.

Return type:

dict[str, NDArray[Any]]

Returns:

dict – Same keys as simulate(), but each value is the concatenation (in ray order, then bounce order) of every ray’s results, shape (n_rays * n_bounces,).

abstractmethod trajectory_segments(pos, vel, n_bounces)[source]#

Trace a trajectory and return its full real-space polyline.

Parameters:

  • pos (ArrayLike) – Initial position (x, y); must lie in the billiard’s interior.

  • vel (ArrayLike) – Initial velocity direction (vx, vy); normalized internally.

  • n_bounces (int) – Number of reflections to simulate.

Return type:

tuple[NDArray[double], dict[str, NDArray[Any]]]

Returns:

  • path (ndarray of float, shape (n_bounces + 1, 2)) – Positions visited, including the initial position as row 0.

  • result (dict) – The same dictionary returned by simulate().

class physicskit.chaos.BunimovichStadium(radius=1.0, straight_length=2.0)[source]#

Bases: _RayTracingBilliard

Chaotic (defocusing) Bunimovich stadium billiard.

Two semicircles joined by straight edges.

Parameters:

  • radius (float) – Radius of the two semicircular end-caps.

  • straight_length (float) – Length of the straight edges joining the semicircles.

Variables:

  • radius (float) – Radius of the two semicircular end-caps.

  • straight_length (float) – Length of the straight edges joining the semicircles.

sample_interior_point()[source]#

Return a point guaranteed to lie in the billiard’s interior.

Return type:

NDArray[double]

Returns:

ndarray of float, shape (2,) – An (x, y) point strictly inside the billiard.

exception physicskit.chaos.ChaoskitError[source]#

Bases: Exception

Base class for all exceptions raised intentionally by physicskit.chaos.

class physicskit.chaos.Chua(alpha=15.6, beta=28.0, m0=-1.1428571428571428, m1=-0.7142857142857143)[source]#

Bases: DynamicalSystem

Chua’s circuit: a simple chaotic electronic oscillator.

Built from just a resistor, two capacitors, an inductor, and one piecewise-linear nonlinear resistor (the “Chua diode”, h below), this is one of the simplest physical systems known to be chaotic, and the first to have its chaos confirmed experimentally in real hardware. For the classic parameters below it produces the famous double-scroll attractor: two spiral lobes, with the trajectory unpredictably switching between them.

\[\begin{split}\dot{x} &= \alpha (y - x - h(x)) \\ \dot{y} &= x - y + z \\ \dot{z} &= -\beta y \\ h(x) &= m_1 x + \tfrac{1}{2}(m_0 - m_1)(|x + 1| - |x - 1|)\end{split}\]
Parameters:

  • alpha (float) – Ratio of the two capacitances.

  • beta (float) – Ratio involving the inductance and second capacitance.

  • m0 (float) – Inner (small-|x|) slope of the Chua diode’s piecewise-linear characteristic.

  • m1 (float) – Outer (large-|x|) slope of the Chua diode’s piecewise-linear characteristic.

Variables:

m1 (alpha, beta, m0,) – System parameters.

dim: int = 3#

State dimension, always 3.

initial_state()[source]#

Default initial condition, slightly off the unstable origin.

Return type:

NDArray[double]

Returns:

ndarray of float, shape (3,)

property params: NDArray[float64]#

Parameter vector (alpha, beta, m0, m1).

Returns:

ndarray of float, shape (4,)

rhs(state, t)[source]#

Evaluate Chua’s circuit vector field.

Parameters:

  • state (NDArray[double]) – State vector (x, y, z).

  • t (float) – Current time (unused; the system is autonomous).

Return type:

NDArray[double]

Returns:

ndarray of float, shape (3,) – Time derivative (dx/dt, dy/dt, dz/dt).

trajectory(state0=None, t0=0.0, dt=0.01, n_steps=20000)[source]#

Integrate a trajectory with the Numba-accelerated RK4 integrator.

Parameters:
Return type:

tuple[NDArray[double], NDArray[double]]

Returns:

  • times (ndarray of float, shape (n_steps + 1,))

  • states (ndarray of float, shape (n_steps + 1, 3))

class physicskit.chaos.CircleBilliard(radius=1.0)[source]#

Bases: _RayTracingBilliard

Integrable circular billiard centered at the origin.

Parameters:

radius (float) – Circle radius.

Variables:

radius (float) – Circle radius.

sample_interior_point()[source]#

Return a point guaranteed to lie in the billiard’s interior.

Return type:

NDArray[double]

Returns:

ndarray of float, shape (2,) – An (x, y) point strictly inside the billiard.

class physicskit.chaos.DiscreteMap[source]#

Bases: ABC

Base class for discrete-time dynamical systems (iterated maps).

dim: int#

Dimension of the state vector x. Set by each concrete subclass.

initial_state()[source]#

Return a reasonable default initial condition, if defined by the subclass.

Return type:

NDArray[double]

Returns:

ndarray of float, shape (dim,) – A default initial state.

Raises:

NotImplementedError – If the subclass does not define a default initial condition.

abstractmethod step(state)[source]#

Advance state by one iteration of the map.

Parameters:

state (NDArray[double]) – Current state vector.

Return type:

NDArray[double]

Returns:

ndarray of float, shape (dim,) – The next state vector.

trajectory(state0=None, n_iter=1000)[source]#

Iterate the map n_iter times starting from state0.

Parameters:
Return type:

NDArray[double]

Returns:

ndarray of float, shape (n_iter + 1, dim) – The state at each iteration, including state0 as row 0.

class physicskit.chaos.DoublePendulum(m1=1.0, m2=1.0, l1=1.0, l2=1.0, g=9.81)[source]#

Bases: DynamicalSystem

Planar double pendulum: point masses on massless rods.

Parameters:

  • m1 (float) – Mass of the first (inner) bob.

  • m2 (float) – Mass of the second (outer) bob.

  • l1 (float) – Length of the first (inner) rod.

  • l2 (float) – Length of the second (outer) rod.

  • g (float) – Gravitational acceleration.

Variables:

g (m1, m2, l1, l2,) – System parameters.

dim: int = 4#

State dimension, always 4. State is (theta1, theta2, omega1, omega2).

energy(state)[source]#

Total mechanical energy, useful for checking integrator drift.

Parameters:

state (NDArray[double]) – State vector (theta1, theta2, omega1, omega2).

Return type:

float

Returns:

float – Total (kinetic + potential) mechanical energy.

initial_state()[source]#

Default initial condition: both rods horizontal, at rest.

Return type:

NDArray[double]

Returns:

ndarray of float, shape (4,)

property params: NDArray[float64]#

Parameter vector (m1, m2, l1, l2, g).

Returns:

ndarray of float, shape (5,)

rhs(state, t)[source]#

Evaluate the double-pendulum vector field.

Parameters:

  • state (NDArray[double]) – State vector (theta1, theta2, omega1, omega2).

  • t (float) – Current time (unused; the system is autonomous).

Return type:

NDArray[double]

Returns:

ndarray of float, shape (4,) – Time derivative (dtheta1/dt, dtheta2/dt, domega1/dt, domega2/dt).

trajectory(state0=None, t0=0.0, dt=0.005, n_steps=10000)[source]#

Integrate a trajectory with the Numba-accelerated RK4 integrator.

Parameters:
Return type:

tuple[NDArray[double], NDArray[double]]

Returns:

  • times (ndarray of float, shape (n_steps + 1,))

  • states (ndarray of float, shape (n_steps + 1, 4))

class physicskit.chaos.Duffing(delta=0.3, alpha=-1.0, beta=1.0, gamma=0.37, omega=1.2)[source]#

Bases: DynamicalSystem

The forced, damped Duffing oscillator.

Governed by x'' + delta*x' + alpha*x + beta*x^3 = gamma*cos(omega*t).

Parameters:

  • delta (float) – Damping coefficient.

  • alpha (float) – Linear stiffness.

  • beta (float) – Cubic (nonlinear) stiffness.

  • gamma (float) – Forcing amplitude.

  • omega (float) – Forcing angular frequency.

Variables:

omega (delta, alpha, beta, gamma,) – System parameters.

dim: int = 2#

State dimension, always 2. State is (x, v).

initial_state()[source]#

Default initial condition (1, 0).

Return type:

NDArray[double]

Returns:

ndarray of float, shape (2,)

property params: NDArray[float64]#

Parameter vector (delta, alpha, beta, gamma, omega).

Returns:

ndarray of float, shape (5,)

rhs(state, t)[source]#

Evaluate the Duffing vector field.

Parameters:

  • state (NDArray[double]) – State vector (x, v).

  • t (float) – Current time (the forcing term depends on t).

Return type:

NDArray[double]

Returns:

ndarray of float, shape (2,) – Time derivative (dx/dt, dv/dt).

trajectory(state0=None, t0=0.0, dt=0.01, n_steps=10000)[source]#

Integrate a trajectory with the Numba-accelerated RK4 integrator.

Parameters:
Return type:

tuple[NDArray[double], NDArray[double]]

Returns:

  • times (ndarray of float, shape (n_steps + 1,))

  • states (ndarray of float, shape (n_steps + 1, 2))

class physicskit.chaos.DynamicalSystem[source]#

Bases: ABC

Base class for continuous-time dynamical systems integrated as ODEs.

Subclasses expose the right-hand side of dx/dt = f(x, t) both as a plain Python method (for convenience/plotting) and, where available, as a Numba-jitted module-level function usable with physicskit.chaos.core.integrators.

dim: int#

Dimension of the state vector x. Set by each concrete subclass.

initial_state()[source]#

Return a reasonable default initial condition, if defined by the subclass.

Return type:

NDArray[double]

Returns:

ndarray of float, shape (dim,) – A default initial state.

Raises:

NotImplementedError – If the subclass does not define a default initial condition.

abstractmethod rhs(state, t)[source]#

Evaluate the vector field at state and time t.

Parameters:
Return type:

NDArray[double]

Returns:

ndarray of float, shape (dim,) – The time derivative dx/dt evaluated at (state, t).

class physicskit.chaos.EllipseBilliard(semi_major=1.5, semi_minor=1.0, n_segments=2000)[source]#

Bases: _RayTracingBilliard

Integrable elliptical billiard, centered at the origin.

The ellipse billiard is integrable: every trajectory remains tangent to a single confocal caustic for all time – either a confocal ellipse (for trajectories that never cross the segment joining the two foci) or a confocal hyperbola (for trajectories that do) – so the Poincare section is foliated by smooth invariant curves, like the Circle and Rectangle billiards. Unlike those two, its boundary is not built from exact circular arcs, so it is represented as a fine closed polygon of n_segments straight edges sampled from the ellipse’s parametric form; the resulting discretization error in the physics is negligible at the default resolution (perimeter error O(1/n_segments^2)).

Parameters:

  • semi_major (float) – Semi-major axis, along x.

  • semi_minor (float) – Semi-minor axis, along y; must be smaller than semi_major.

  • n_segments (int) – Number of straight edges used to approximate the smooth ellipse.

Variables:

  • semi_minor (semi_major,) – Ellipse semi-axes.

  • n_segments (int) – Polygon-approximation resolution.

Raises:

ValueError – If semi_major or semi_minor is not positive, or semi_major does not exceed semi_minor.

foci()[source]#

The two foci of the ellipse, at (+-c, 0) with c = sqrt(a^2 - b^2).

Return type:

tuple[NDArray[double], NDArray[double]]

Returns:

f1, f2 (ndarray of float, shape (2,)) – The two focal points.

sample_interior_point()[source]#

Return a point guaranteed to lie in the billiard’s interior.

Return type:

NDArray[double]

Returns:

ndarray of float, shape (2,) – An (x, y) point strictly inside the billiard.

class physicskit.chaos.HenonMap(a=1.4, b=0.3)[source]#

Bases: DiscreteMap

The Henon map x' = 1 - a*x^2 + y, y' = b*x.

The classic chaotic parameters are a=1.4, b=0.3.

Parameters:

  • a (float) – Map parameter.

  • b (float) – Map parameter.

Variables:

b (a,) – Map parameters.

dim: int = 2#

State dimension, always 2. State is (x, y).

initial_state()[source]#

Default initial condition (0, 0).

Return type:

NDArray[double]

Returns:

ndarray of float, shape (2,)

step(state)[source]#

Advance (x, y) by one map iteration.

Parameters:

state (NDArray[double]) – Current state (x, y).

Return type:

NDArray[double]

Returns:

ndarray of float, shape (2,) – Next state (x, y).

trajectory(state0=None, n_iter=1000)[source]#

Iterate the map n_iter times starting from state0.

Parameters:
Return type:

NDArray[double]

Returns:

ndarray of float, shape (n_iter + 1, 2) – State at each iteration, including state0 as row 0.

exception physicskit.chaos.InvalidParameterError[source]#

Bases: ChaoskitError, ValueError

A system or tool was constructed or called with an invalid parameter.

Examples include a billiard’s scatterer not fitting inside its cell, or a map’s parameter falling outside its valid range.

class physicskit.chaos.LogisticMap(r=3.9)[source]#

Bases: DiscreteMap

The logistic map x' = r*x*(1-x): the simplest gateway to chaos.

Varying the single growth-rate parameter r takes this map through the complete period-doubling route to chaos: a stable fixed point for r < 3, then successive period-doubling bifurcations at r values that accumulate geometrically (ratio converging to the universal Feigenbaum constant delta ~= 4.669) onto the onset of chaos at r ~= 3.5699, beyond which the map is chaotic for most (but not all – note the periodic windows, the largest at r ~= 3.8284) values of r up to r=4.

Parameters:

r (float) – Growth rate parameter; interesting values range over [0, 4].

Variables:

r (float) – Growth rate parameter.

dim: int = 1#

State dimension, always 1. State is (x,).

initial_state()[source]#

Default initial condition (0.5,).

Return type:

NDArray[double]

Returns:

ndarray of float, shape (1,)

step(state)[source]#

Advance x by one map iteration.

Parameters:

state (NDArray[double]) – Current state (x,).

Return type:

NDArray[double]

Returns:

ndarray of float, shape (1,) – Next state (x,).

trajectory(state0=None, n_iter=1000)[source]#

Iterate the map n_iter times starting from state0.

Parameters:
Return type:

NDArray[double]

Returns:

ndarray of float, shape (n_iter + 1, 1) – State at each iteration, including state0 as row 0.

class physicskit.chaos.Lorenz(sigma=10.0, rho=28.0, beta=2.6666666666666665)[source]#

Bases: DynamicalSystem

The Lorenz attractor.

Parameters:

  • sigma (float) – Prandtl-number-like parameter.

  • rho (float) – Rayleigh-number-like parameter.

  • beta (float) – Geometric parameter.

Variables:

beta (sigma, rho,) – System parameters.

dim: int = 3#

State dimension, always 3.

initial_state()[source]#

Default initial condition (1, 1, 1).

Return type:

NDArray[double]

Returns:

ndarray of float, shape (3,)

property params: NDArray[float64]#

Parameter vector (sigma, rho, beta).

Returns:

ndarray of float, shape (3,)

rhs(state, t)[source]#

Evaluate the Lorenz vector field.

Parameters:

  • state (NDArray[double]) – State vector (x, y, z).

  • t (float) – Current time (unused; the system is autonomous).

Return type:

NDArray[double]

Returns:

ndarray of float, shape (3,) – Time derivative (dx/dt, dy/dt, dz/dt).

trajectory(state0=None, t0=0.0, dt=0.01, n_steps=10000)[source]#

Integrate a trajectory with the Numba-accelerated RK4 integrator.

Parameters:
Return type:

tuple[NDArray[double], NDArray[double]]

Returns:

  • times (ndarray of float, shape (n_steps + 1,))

  • states (ndarray of float, shape (n_steps + 1, 3))

class physicskit.chaos.MagneticPendulum(magnet_positions=None, friction=0.2, spring=0.2, height=0.2, strength=1.0)[source]#

Bases: DynamicalSystem

A pendulum bob swinging over several fixed magnets: a multistable, chaotic system.

A damped pendulum bob, modeled in the small-swing (flat, 2D) limit, is pulled down toward the origin by a linear restoring force and attracted toward each of several fixed magnets by an inverse-square-like force (softened by a “height” offset height, the bob’s height above the magnet plane, which avoids a force singularity directly above a magnet). Friction eventually settles the bob at rest near whichever magnet “won” – but which magnet wins depends on the starting position with famously fractal sensitivity, making this the classic system for visualizing fractal basin boundaries; see physicskit.chaos.visualizers.basins.plot_basin_of_attraction().

Parameters:

  • magnet_positions (NDArray[double] | None) – Magnet (x, y) positions; defaults to 3 magnets at the vertices of an equilateral triangle inscribed in the unit circle.

  • friction (float) – Damping coefficient.

  • spring (float) – Linear restoring-force coefficient (pulling the bob back toward the origin, as in the small-swing limit of gravity).

  • height (float) – The bob’s height above the magnet plane; softens the force near a magnet (larger values give a gentler, less singular pull).

  • strength (float) – Magnet attraction strength.

Variables:

  • magnet_positions (ndarray of float, shape (n_magnets, 2)) – Magnet positions.

  • strength (friction, spring, height,) – System parameters.

dim: int = 4#

State dimension, always 4. State is (x, y, vx, vy).

initial_state()[source]#

Default initial condition: released from rest, off-center.

Return type:

NDArray[double]

Returns:

ndarray of float, shape (4,)

property params: NDArray[float64]#

Parameter vector (friction, spring, height, strength, mx_0, my_0, ...).

Returns:

ndarray of float, shape (4 + 2*n_magnets,)

rhs(state, t)[source]#

Evaluate the magnetic pendulum vector field.

Parameters:

  • state (NDArray[double]) – State vector (x, y, vx, vy).

  • t (float) – Current time (unused; the system is autonomous).

Return type:

NDArray[double]

Returns:

ndarray of float, shape (4,) – Time derivative (dx/dt, dy/dt, dvx/dt, dvy/dt).

trajectory(state0=None, t0=0.0, dt=0.02, n_steps=5000)[source]#

Integrate a trajectory with the Numba-accelerated RK4 integrator.

Parameters:
Return type:

tuple[NDArray[double], NDArray[double]]

Returns:

  • times (ndarray of float, shape (n_steps + 1,))

  • states (ndarray of float, shape (n_steps + 1, 4))

class physicskit.chaos.QuantumBakersMap(dim, alpha=0.5)[source]#

Bases: object

The quantum baker’s map: the (generalized) baker’s map’s quantization.

Built via the Balazs-Voros/Saraceno construction, generalized to an arbitrary cut alpha in exact correspondence with physicskit.chaos’s classical BakersMap: the Floquet operator applies the discrete Fourier transform separately to the q < alpha and q >= alpha position-basis blocks (mirroring the classical map stretching each piece independently), then transforms the result back to the full position representation.

Parameters:

  • dim (int) – Hilbert space dimension (number of position basis states).

  • alpha (float) – Cut position, in (0, 1); matches the classical BakersMap’s alpha. The classic Balazs-Voros construction is the alpha=0.5 case with dim even.

Variables:

  • dim (int) – Hilbert space dimension.

  • alpha (float) – Cut position.

  • hbar (float) – Effective reduced Planck constant, 1 / (2*pi*dim): the unit torus holds dim = 1/h = 1/(2*pi*hbar) states (Balazs & Voros 1989), so plane waves are exp(i*p*q/hbar) = exp(2*pi*i*dim*p*q).

Raises:

InvalidParameterError – If alpha does not satisfy 0 < alpha < 1, or if rounding alpha * dim to the nearest integer would leave either of the two blocks with fewer than 1 basis state.

Notes

Exact correspondence with the classical cut requires alpha * dim to be an integer; for other values, the nearest integer split is used, and alpha is left at the value the caller requested rather than silently adjusted to the value actually realized – pass a dim that makes alpha * dim (near-)integral for the closest match.

coherent_state(q0, p0)[source]#

A minimum-uncertainty wavepacket centered at (q0, p0).

Parameters:

  • q0 (float) – Phase-space center of the wavepacket, each in [0, 1).

  • p0 (float) – Phase-space center of the wavepacket, each in [0, 1).

Return type:

NDArray[cdouble]

Returns:

ndarray of complex, shape (dim,) – Normalized coherent state in the position representation.

eigenphases()[source]#

Quasi-energies (eigenphases) of the Floquet operator.

Hand these to physicskit.rmt to study their spacing statistics.

Return type:

NDArray[double]

Returns:

ndarray of float, shape (dim,) – Eigenphases, sorted ascending, in radians in (-pi, pi].

evolve(psi, n_steps=1)[source]#

Propagate a state through n_steps map iterations.

Parameters:

  • psi (ArrayLike) – Initial state in the position representation; normalized internally.

  • n_steps (int) – Number of map iterations to advance.

Return type:

NDArray[cdouble]

Returns:

ndarray of complex, shape (n_steps + 1, dim) – The state after each iteration, including the (normalized) initial state as row 0.

floquet_operator()[source]#

Build the one-iteration Floquet (evolution) operator.

Return type:

NDArray[cdouble]

Returns:

ndarray of complex, shape (dim, dim) – Unitary matrix advancing a state, in the position representation, by one map iteration.

husimi(psi, resolution=80, n_wraps=3)[source]#

Husimi phase-space distribution of a state; see husimi_function().

Parameters:

  • psi (ArrayLike) – State in the position representation.

  • resolution (int) – Number of grid points along each of the q and p axes.

  • n_wraps (int) – Number of periodic images summed to periodize the coherent states used.

Return type:

tuple[NDArray[double], NDArray[double], NDArray[double]]

Returns:

Q, P, husimi (ndarray of float, shape (resolution, resolution)) – Phase-space grid and the (peak-normalized) Husimi distribution.

class physicskit.chaos.QuantumBilliard(billiard, resolution=150)[source]#

Bases: object

Dirichlet Helmholtz eigenstates of a billiard: “particle in a box” quantum chaos.

A quantum particle confined to a chaotic billiard is one of the two textbook playgrounds of quantum chaos (alongside quantized maps like QuantumKickedRotor): its energy eigenvalues E_n = k_n^2 (in units where hbar^2 / 2m = 1) obey Weyl’s law on average (see weyl_counting_function()) but fluctuate around it in a way that reflects the underlying classical dynamics, and its eigenfunctions can “scar” – show anomalously enhanced density – on unstable classical periodic orbits.

Built directly from any BilliardSystem already in physicskit.chaos (no shape-specific code needed): a regular grid is laid over the shape’s bounding box, points inside the boundary (via points_in_billiard()) become unknowns of a standard five-point finite-difference Laplacian, and the lowest eigenpairs of the resulting sparse, symmetric positive-definite matrix are found by shift-invert Lanczos iteration.

Parameters:

  • billiard (BilliardSystem) – The billiard shape to quantize.

  • resolution (int) – Number of grid points along the longer side of the billiard’s bounding box; the grid spacing (and hence both the accuracy and the cost of solving for eigenstates) scales with this.

Variables:

  • billiard (BilliardSystem) – The billiard shape.

  • resolution (int) – Grid resolution.

Raises:

InvalidParameterError – If resolution is smaller than 10.

Notes

Finite differences converge slowly (error O(h^2)) and, being defined on a Cartesian grid, represent curved or slanted boundaries only approximately; treat eigenvalues as accurate to a few percent at the default resolution; increase resolution for tighter results, at roughly quadratic cost in memory and eigensolver runtime.

area()[source]#

Interior area, estimated by counting grid points inside the boundary.

Return type:

float

Returns:

float – Approximate billiard area.

eigenstates(n_states=6)[source]#

Solve for the lowest n_states Dirichlet eigenpairs.

Parameters:

n_states (int) – Number of lowest eigenstates to compute.

Return type:

tuple[NDArray[double], NDArray[double]]

Returns:

  • eigenvalues (ndarray of float, shape (n_states,)) – Eigenvalues k_n^2, ascending.

  • eigenfunctions (ndarray of float, shape (n_states, nx, ny)) – Eigenfunctions on the grid returned by grid(), each normalized to a peak absolute value of 1; grid points outside the billiard are nan (so plt.imshow/pcolormesh leave them blank).

grid()[source]#

The Cartesian grid eigenfunctions are returned on.

Return type:

tuple[NDArray[double], NDArray[double]]

Returns:

X, Y (ndarray of float, shape (nx, ny)) – Grid coordinates (as from np.meshgrid(..., indexing="ij")).

wavenumbers(n_states=6)[source]#

Wavenumbers k_n = sqrt(eigenvalue) of the lowest n_states states.

Parameters:

n_states (int) – Number of lowest eigenstates to compute.

Return type:

NDArray[double]

Returns:

ndarray of float, shape (n_states,) – Wavenumbers, ascending.

weyl_counting_function(k)[source]#

Weyl’s law: the average number of eigenvalues below wavenumber k.

N(k) ~ Area * k^2 / (4*pi) - Perimeter * k / (4*pi), the leading area term plus the first-order (Dirichlet) boundary correction. Real billiards’ actual eigenvalue counts fluctuate around this smooth curve; how they fluctuate (level “rigidity”) is a hallmark of whether the classical billiard is integrable or chaotic – exactly the kind of spectral statistic physicskit.rmt is built to quantify, using wavenumbers() as its input.

Parameters:

k (NDArray[double] | float) – Wavenumber(s) at which to evaluate the counting function.

Return type:

NDArray[double] | float

Returns:

float or ndarray of float – Expected (smoothed) number of eigenvalues below k.

class physicskit.chaos.QuantumKickedRotor(k=1.0, dim=64, hbar=None)[source]#

Bases: object

The quantum kicked rotor: the standard map’s quantization.

Built as the one-period Floquet operator of a rotor periodically kicked by a potential k * cos(theta), in exact correspondence with physicskit.chaos’s classical StandardMap (p_new = p + k*sin(theta), theta_new = theta + p_new): the Hilbert space is the dim-point position (angle) representation on [0, 2*pi), and the Floquet operator alternates a kick phase (diagonal in the angle basis) with a free-rotation phase (diagonal in the momentum basis), transforming between the two via the discrete Fourier transform.

As dim grows (equivalently, as hbar shrinks towards its default 2*pi/dim), the quantum dynamics of a narrow wavepacket increasingly tracks the corresponding classical StandardMap orbit, until the packet spreads across a chaotic region – the quantum-classical correspondence breaking down being one of the central phenomena of quantum chaos.

Parameters:

  • k (float) – Kick strength; matches the classical StandardMap’s k exactly.

  • dim (int) – Hilbert space dimension (number of angle basis states).

  • hbar (float | None) – Effective Planck constant; defaults to 2*pi/dim, the standard choice that keeps the quantized torus’s phase-space cell count equal to dim.

Variables:

  • k (float) – Kick strength.

  • dim (int) – Hilbert space dimension.

  • hbar (float) – Effective Planck constant.

Raises:

InvalidParameterError – If dim is smaller than 2.

coherent_state(theta0, p0)[source]#

A minimum-uncertainty wavepacket centered at (theta0, p0).

Useful as a semiclassical initial state for evolve(), to watch the quantum dynamics track (and eventually depart from) the corresponding classical orbit.

Parameters:

  • theta0 (float) – Phase-space center of the wavepacket.

  • p0 (float) – Phase-space center of the wavepacket.

Return type:

NDArray[cdouble]

Returns:

ndarray of complex, shape (dim,) – Normalized coherent state in the angle representation.

eigenphases()[source]#

Quasi-energies (eigenphases) of the Floquet operator.

Hand these to physicskit.rmt to study their spacing statistics.

Return type:

NDArray[double]

Returns:

ndarray of float, shape (dim,) – Eigenphases (angles of the Floquet operator’s unit-modulus eigenvalues), sorted ascending, in radians in (-pi, pi].

evolve(psi, n_steps=1)[source]#

Propagate a state through n_steps kicks.

Parameters:

  • psi (ArrayLike) – Initial state in the angle representation; normalized internally.

  • n_steps (int) – Number of Floquet periods to advance.

Return type:

NDArray[cdouble]

Returns:

ndarray of complex, shape (n_steps + 1, dim) – The state after each kick, including the (normalized) initial state as row 0.

floquet_operator()[source]#

Build the one-period Floquet (evolution) operator.

Return type:

NDArray[cdouble]

Returns:

ndarray of complex, shape (dim, dim) – Unitary matrix advancing a state, in the angle representation, by one kick-and-rotation period.

husimi(psi, resolution=80, n_wraps=3)[source]#

Husimi phase-space distribution of a state; see husimi_function().

Parameters:

  • psi (ArrayLike) – State in the angle representation.

  • resolution (int) – Number of grid points along each of the theta and p axes.

  • n_wraps (int) – Number of periodic images summed to periodize the coherent states used.

Return type:

tuple[NDArray[double], NDArray[double], NDArray[double]]

Returns:

Theta, P, husimi (ndarray of float, shape (resolution, resolution)) – Phase-space grid and the (peak-normalized) Husimi distribution.

class physicskit.chaos.RectangleBilliard(width=2.0, height=1.0)[source]#

Bases: _RayTracingBilliard

Integrable rectangular billiard centered at the origin.

Parameters:

  • width (float) – Full extent along x.

  • height (float) – Full extent along y.

Variables:

  • width (float) – Full extent along x.

  • height (float) – Full extent along y.

sample_interior_point()[source]#

Return a point guaranteed to lie in the billiard’s interior.

Return type:

NDArray[double]

Returns:

ndarray of float, shape (2,) – An (x, y) point strictly inside the billiard.

class physicskit.chaos.RestrictedThreeBody(mu=0.012277471)[source]#

Bases: DynamicalSystem

The planar circular restricted three-body problem (CR3BP).

A massless third body moves under the gravity of two massive primaries (masses 1 - mu and mu, in normalized units) that are themselves in a fixed circular orbit about their common center of mass. In the rotating (co-precessing) reference frame, the primaries sit fixed at (-mu, 0) and (1 - mu, 0), and the third body’s motion picks up centrifugal and Coriolis terms alongside the two gravitational pulls.

This system is what led Poincare to the first discovery of deterministic chaos: he found that, unlike the exactly solvable two-body problem, CR3BP trajectories can depend on initial conditions in an essentially unpredictable way. It also genuinely coexists with regular (quasi-periodic, KAM-stable) motion – the default initial condition below is the classic Arenstorf orbit, a stable periodic orbit famous in the numerical-methods literature as an ODE-solver stress test (it passes very close to the smaller primary); see the example gallery for a nearby, only slightly perturbed initial condition that is chaotic instead.

Because of the velocity-dependent Coriolis terms, this system is not a separable Hamiltonian of the form pos'' = force(pos, t), so it must be integrated with physicskit.chaos.core.integrators.rk4_integrate() (the symplectic leapfrog_integrate() / yoshida4_integrate() do not apply).

Parameters:

mu (float) – Mass parameter (mass of the smaller primary, in units where the total mass is 1); the default is the Earth-Moon-like value used in the classic Arenstorf orbit.

Variables:

mu (float) – Mass parameter.

dim: int = 4#

State dimension, always 4. State is (x, y, vx, vy).

initial_state()[source]#

The classic Arenstorf periodic-orbit initial condition.

Return type:

NDArray[double]

Returns:

ndarray of float, shape (4,)

jacobi_constant(state)[source]#

The Jacobi constant: the CR3BP’s conserved rotating-frame energy analog.

Parameters:

state (NDArray[double]) – State vector (x, y, vx, vy).

Return type:

float

Returns:

float – The Jacobi constant C = 2*Omega(x, y) - (vx^2 + vy^2), where Omega is the effective (gravitational + centrifugal) potential. Conserved along any trajectory; useful for checking integrator fidelity (e.g. with physicskit.chaos.utils.metrics.energy_drift()).

property params: NDArray[float64]#

Parameter vector (mu,).

Returns:

ndarray of float, shape (1,)

rhs(state, t)[source]#

Evaluate the CR3BP vector field.

Parameters:

  • state (NDArray[double]) – State vector (x, y, vx, vy).

  • t (float) – Current time (unused; the system is autonomous in the rotating frame).

Return type:

NDArray[double]

Returns:

ndarray of float, shape (4,) – Time derivative (dx/dt, dy/dt, dvx/dt, dvy/dt).

trajectory(state0=None, t0=0.0, dt=0.0005, n_steps=40000)[source]#

Integrate a trajectory with the Numba-accelerated RK4 integrator.

Parameters:

  • state0 (NDArray[double] | None) – Initial state; defaults to initial_state().

  • t0 (float) – Initial time.

  • dt (float) – Integration step size (small, since the Arenstorf orbit’s default initial condition passes very close to the smaller primary).

  • n_steps (int) – Number of integration steps.

Return type:

tuple[NDArray[double], NDArray[double]]

Returns:

  • times (ndarray of float, shape (n_steps + 1,))

  • states (ndarray of float, shape (n_steps + 1, 4))

class physicskit.chaos.Rossler(a=0.2, b=0.2, c=5.7)[source]#

Bases: DynamicalSystem

The Rossler attractor.

Parameters:

  • a (float) – System parameter.

  • b (float) – System parameter.

  • c (float) – System parameter.

Variables:

c (a, b,) – System parameters.

dim: int = 3#

State dimension, always 3.

initial_state()[source]#

Default initial condition (1, 1, 1).

Return type:

NDArray[double]

Returns:

ndarray of float, shape (3,)

property params: NDArray[float64]#

Parameter vector (a, b, c).

Returns:

ndarray of float, shape (3,)

rhs(state, t)[source]#

Evaluate the Rossler vector field.

Parameters:

  • state (NDArray[double]) – State vector (x, y, z).

  • t (float) – Current time (unused; the system is autonomous).

Return type:

NDArray[double]

Returns:

ndarray of float, shape (3,) – Time derivative (dx/dt, dy/dt, dz/dt).

trajectory(state0=None, t0=0.0, dt=0.01, n_steps=10000)[source]#

Integrate a trajectory with the Numba-accelerated RK4 integrator.

Parameters:
Return type:

tuple[NDArray[double], NDArray[double]]

Returns:

  • times (ndarray of float, shape (n_steps + 1,))

  • states (ndarray of float, shape (n_steps + 1, 3))

class physicskit.chaos.SinaiBilliard(cell_size=2.0, scatterer_radius=0.5)[source]#

Bases: _RayTracingBilliard

Chaotic (defocusing) Sinai billiard.

A square cell with a circular scatterer removed from its center.

Parameters:

  • cell_size (float) – Full side length of the square cell.

  • scatterer_radius (float) – Radius of the central circular scatterer; must be smaller than half of cell_size.

Variables:

  • cell_size (float) – Full side length of the square cell.

  • scatterer_radius (float) – Radius of the central circular scatterer.

Raises:

ValueError – If scatterer_radius is not smaller than half of cell_size.

sample_interior_point()[source]#

Return a point guaranteed to lie in the billiard’s interior.

Return type:

NDArray[double]

Returns:

ndarray of float, shape (2,) – An (x, y) point strictly inside the billiard.

class physicskit.chaos.StandardMap(k=1.0)[source]#

Bases: DiscreteMap

The Chirikov-Taylor standard map on the (theta, p) cylinder.

k=0 is integrable; chaos onset is around k~1, with global chaos for k >~ 4-5.

Parameters:

k (float) – Kick strength.

Variables:

k (float) – Kick strength.

dim: int = 2#

State dimension, always 2. State is (theta, p).

initial_state()[source]#

Default initial condition (0.1, 0.1).

Return type:

NDArray[double]

Returns:

ndarray of float, shape (2,)

step(state)[source]#

Advance (theta, p) by one map iteration.

Parameters:

state (NDArray[double]) – Current state (theta, p).

Return type:

NDArray[double]

Returns:

ndarray of float, shape (2,) – Next state (theta, p).

trajectory(state0=None, n_iter=1000)[source]#

Iterate the map n_iter times starting from state0.

Parameters:
Return type:

NDArray[double]

Returns:

ndarray of float, shape (n_iter + 1, 2) – State at each iteration, including state0 as row 0.

class physicskit.chaos.TruncatedCircleBilliard(radius=1.0, cut=0.3)[source]#

Bases: _RayTracingBilliard

A disk truncated by a straight chord.

The region x <= radius - cut of the disk of the given radius is kept.

Parameters:

  • radius (float) – Disk radius.

  • cut (float) – Distance the chord is cut in from the disk’s edge; must satisfy 0 < cut < radius.

Variables:

  • radius (float) – Disk radius.

  • cut (float) – Distance the chord is cut in from the disk’s edge.

Raises:

ValueError – If cut does not satisfy 0 < cut < radius.

sample_interior_point()[source]#

Return a point guaranteed to lie in the billiard’s interior.

Return type:

NDArray[double]

Returns:

ndarray of float, shape (2,) – An (x, y) point strictly inside the billiard.

On this page
Show Source