physicskit.chaos#
physicskit.chaos: visual analysis and simulation of chaotic dynamical systems and 2D billiards.
- class physicskit.chaos.BakersMap(alpha=0.5)[source]#
Bases:
DiscreteMapThe (generalized) baker’s map on the unit square: the textbook model of chaos.
The map cuts the unit square
[0, 1) x [0, 1)atx = 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 hasalpha=0.5.- Variables:
alpha (float) – Cut position.
- Raises:
ValueError – If alpha does not satisfy
0 < alpha < 1.
Notes
x = 0is an exact fixed point of the map for every alpha (since0 / 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 exactly0.0in 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 ontrajectory()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;HenonMapand continuous systems such asLorenzdo not share this failure mode and are better suited to long-trajectory statistics.- 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/alphabranch and1 - alphain the expanding-by1/(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 validatingphysicskit.chaos.utils.metrics.lyapunov_exponent_from_divergence()and similar numerical estimators.
- class physicskit.chaos.BilliardSystem[source]#
Bases:
ABCBase 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 (arclengths, sine of the reflection anglesin(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:
- 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:
- 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:
- Returns:
float – The billiard boundary’s total arclength.
- abstractmethod sample_interior_point()[source]#
Return a point guaranteed to lie in the billiard’s interior.
- abstractmethod simulate(pos, vel, n_bounces)[source]#
Trace
n_bouncesspecular reflections from an initial ray.- Parameters:
- Return type:
- 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:
- Return type:
- 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:
- Return type:
- 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.Brusselator(a=1.0, b=2.5)[source]#
Bases:
DynamicalSystemThe Brusselator chemical oscillator of Prigogine and Lefever (1968).
\[\dot x = a - (b + 1)x + x^2 y, \qquad \dot y = bx - x^2 y.\]The unique fixed point \((a, b/a)\) has Jacobian trace \(b - 1 - a^2\) and determinant \(a^2\), so it undergoes a supercritical Hopf bifurcation at \(b_c = 1 + a^2\), with frequency \(\omega = a\) at onset.
Examples
>>> Brusselator(a=1.0).hopf_point 2.0
- eigenvalues()[source]#
Eigenvalues of
jacobian().
- property fixed_point: NDArray[float64]#
The fixed point \((a, b/a)\).
- jacobian()[source]#
Jacobian at the fixed point, \(\begin{pmatrix} b - 1 & a^2 \\ -b & -a^2 \end{pmatrix}\).
- property params: NDArray[float64]#
Parameter vector
(a, b).
- class physicskit.chaos.BunimovichStadium(radius=1.0, straight_length=2.0)[source]#
Bases:
_RayTracingBilliardChaotic (defocusing) Bunimovich stadium billiard.
Two semicircles joined by straight edges.
- Parameters:
- Variables:
- exception physicskit.chaos.ChaoskitError[source]#
Bases:
ExceptionBase 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:
DynamicalSystemChua’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”,
hbelow), 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:
- Variables:
- property params: NDArray[float64]#
Parameter vector
(alpha, beta, m0, m1).- Returns:
ndarray of float, shape (4,)
- class physicskit.chaos.CircleBilliard(radius=1.0)[source]#
Bases:
_RayTracingBilliardIntegrable circular billiard centered at the origin.
- class physicskit.chaos.DiscreteMap[source]#
Bases:
ABCBase class for discrete-time dynamical systems (iterated maps).
- initial_state()[source]#
Return a reasonable default initial condition, if defined by the subclass.
- Return type:
- Returns:
ndarray of float, shape (dim,) – A default initial state.
- Raises:
NotImplementedError – If the subclass does not define a default initial condition.
- class physicskit.chaos.DoublePendulum(m1=1.0, m2=1.0, l1=1.0, l2=1.0, g=9.81)[source]#
Bases:
DynamicalSystemPlanar double pendulum: point masses on massless rods.
- Parameters:
- Variables:
g (m1, m2, l1, l2,) – System parameters.
- property params: NDArray[float64]#
Parameter vector
(m1, m2, l1, l2, g).- Returns:
ndarray of float, shape (5,)
- class physicskit.chaos.Duffing(delta=0.3, alpha=-1.0, beta=1.0, gamma=0.37, omega=1.2)[source]#
Bases:
DynamicalSystemThe forced, damped Duffing oscillator.
Governed by
x'' + delta*x' + alpha*x + beta*x^3 = gamma*cos(omega*t).- Parameters:
- Variables:
- property params: NDArray[float64]#
Parameter vector
(delta, alpha, beta, gamma, omega).- Returns:
ndarray of float, shape (5,)
- class physicskit.chaos.DynamicalSystem[source]#
Bases:
ABCBase 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 withphysicskit.chaos.core.integrators.- initial_state()[source]#
Return a reasonable default initial condition, if defined by the subclass.
- Return type:
- Returns:
ndarray of float, shape (dim,) – A default initial state.
- Raises:
NotImplementedError – If the subclass does not define a default initial condition.
- class physicskit.chaos.EllipseBilliard(semi_major=1.5, semi_minor=1.0, n_segments=2000)[source]#
Bases:
_RayTracingBilliardIntegrable 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:
- 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.
- class physicskit.chaos.ForcedVanDerPol(mu=10.0, A=14.5, omega=2.5)[source]#
Bases:
DynamicalSystemThe sinusoidally forced Van der Pol oscillator.
Governed by
x'' - mu*(1 - x^2)*x' + x = A*cos(omega*t): a self-excited (negatively damped at small amplitude) valve-circuit oscillator driven by a periodic signal. WithA = 0it is the plain Van der Pol oscillator, whose limit cycle has amplitude close to 2 for small mu and becomes a relaxation oscillation of period roughly(3 - 2 ln 2) mufor large mu.Cartwright & Littlewood (1945) proved, in the scaling
A = b * omega * muwith mu large, that for a range of b the forced oscillator has two coexisting stable periodic motions whose periods are different odd multiples of the forcing period2*pi/omega, plus an invariant “bad” set of infinitely many periodic and uncountably many non-periodic orbits – the first proof of chaos in an equation from physics. The defaults below are one such parameter set (b = 0.58), where stable subharmonics of period3Tand5Tcoexist.- Parameters:
- Variables:
omega (mu, A,) – System parameters.
References
M. L. Cartwright and J. E. Littlewood, “On Non-Linear Differential Equations of the Second Order: I. The Equation y’’ - k(1 - y^2)y’ + y = b lambda k cos(lambda t + a), k Large,” J. London Math. Soc. 20, 180-189 (1945).
Examples
>>> import numpy as np >>> from physicskit.chaos.systems.continuous import ForcedVanDerPol >>> system = ForcedVanDerPol(mu=1.0, A=0.5, omega=2.0) >>> system.rhs(np.array([0.0, 1.0]), 0.0).tolist() [1.0, 1.5] >>> system.forcing_period == 2 * np.pi / 2.0 True
- property params: NDArray[float64]#
Parameter vector
(mu, A, omega).- Returns:
ndarray of float, shape (3,)
- stroboscopic_map(states0, n_periods, steps_per_period=400)[source]#
Sample many trajectories once per forcing period, in parallel.
Every trajectory starts at
t = 0and is integrated with RK4 at a fixed stepforcing_period / steps_per_period, recording the state att = 0, T, 2T, .... This stroboscopic (Poincare) map turns the periodic orbits of the flow into finite cycles of points – a subharmonic of periodnTbecomes an n-cycle – which is how Cartwright and Littlewood’s coexisting periodic motions are told apart. Many initial conditions (e.g. a basin-of-attraction grid) are integrated at once on all cores.- Parameters:
- Return type:
- Returns:
ndarray of float, shape (n, n_periods + 1, 2) –
out[i, k]is the state of trajectory i att = k*T(out[i, 0]is its initial state). A single(2,)input still returns a leading axis of length 1.
Examples
>>> import numpy as np >>> from physicskit.chaos.systems.continuous import ForcedVanDerPol >>> linear = ForcedVanDerPol(mu=0.0, A=0.0, omega=1.0) >>> out = linear.stroboscopic_map([1.0, 0.0], n_periods=3) >>> out.shape (1, 4, 2) >>> bool(np.allclose(out[0, -1], [1.0, 0.0], atol=1e-8)) # x = cos(t) has period T True
- class physicskit.chaos.HenonMap(a=1.4, b=0.3)[source]#
Bases:
DiscreteMapThe Henon map
x' = 1 - a*x^2 + y, y' = b*x.The classic chaotic parameters are
a=1.4, b=0.3.- Parameters:
- Variables:
b (a,) – Map parameters.
- class physicskit.chaos.HopfNormalForm(mu=0.1, omega=1.0, a=-1.0, c=0.0)[source]#
Bases:
DynamicalSystemNormal form of the Hopf bifurcation, \(\dot r = \mu r + a r^3 + c r^5\), \(\dot\phi = \omega\).
- Parameters:
mu (
float) – Bifurcation parameter: the real part of the fixed point’s eigenvalues.omega (
float) – Rotation frequency at onset.a (
float) – Cubic (first Lyapunov) coefficient; negative is supercritical, positive is subcritical.c (
float) – Quintic coefficient; must be negative whena > 0to keep the flow bounded.
Examples
>>> HopfNormalForm(mu=0.25).limit_cycle_radii() [(0.5, 'stable')] >>> [(round(r, 3), s) for r, s in HopfNormalForm(mu=-0.1, a=1.0, c=-1.0).limit_cycle_radii()] [(0.336, 'unstable'), (0.942, 'stable')]
- property eigenvalues: NDArray[complex128]#
Eigenvalues \(\mu \pm i\omega\) of the fixed point at the origin.
- limit_cycle_radii()[source]#
Radii of the limit cycles and their stability.
The cycles are the positive roots \(r = \sqrt{\rho}\) of \(\mu + a\rho + c\rho^2 = 0\); a cycle is stable when \(\dot r\) decreases through it.
- property params: NDArray[float64]#
Parameter vector
(mu, omega, a, c).
- exception physicskit.chaos.InvalidParameterError[source]#
Bases:
ChaoskitError,ValueErrorA 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.Kuramoto(n_oscillators=500, K=1.0, gamma=0.5, distribution='lorentzian')[source]#
Bases:
DynamicalSystemThe Kuramoto model of \(N\) all-to-all coupled phase oscillators.
Natural frequencies are placed deterministically at the quantiles of \(g(\omega)\), \(\omega_i = G^{-1}((i - 1/2)/N)\), which removes the sample-to-sample scatter of a random draw and makes finite-\(N\) runs follow the \(N \to \infty\) theory closely.
- Parameters:
- Variables:
omega (ndarray of float, shape (n_oscillators,)) – Natural frequencies.
Examples
>>> model = Kuramoto(n_oscillators=400, K=3.0, gamma=0.5) >>> model.critical_coupling 1.0 >>> r = model.order_parameter_series(t_max=60.0, seed=0) >>> bool(abs(r[-200:].mean() - kuramoto_order_parameter_lorentzian(3.0, 0.5)) < 0.03) True
- property critical_coupling: float#
Kuramoto’s onset of synchronization, \(K_c = 2/(\pi g(0))\).
\(2\gamma\) for a Lorentzian and \(\sqrt{8/\pi}\,\sigma\) for a Gaussian of standard deviation \(\sigma\).
- static order_parameter(states)[source]#
Order parameter \(r = |\langle e^{i\theta_j} \rangle_j|\) of each row of phases.
- order_parameter_series(t_max=100.0, dt=0.05, state0=None, seed=0, record_every=1)[source]#
Integrate and return only \(r(t)\), which keeps memory flat in \(N\).
- Parameters:
- Return type:
- Returns:
ndarray of float – \(r\) at times
0, record_every*dt, ....
- property params: NDArray[float64]#
Parameter vector
(K, omega_1, ..., omega_N).
- class physicskit.chaos.LogisticMap(r=3.9)[source]#
Bases:
DiscreteMapThe 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 atrvalues that accumulate geometrically (ratio converging to the universal Feigenbaum constantdelta ~= 4.669) onto the onset of chaos atr ~= 3.5699, beyond which the map is chaotic for most (but not all – note the periodic windows, the largest atr ~= 3.8284) values of r up tor=4.- Parameters:
r (
float) – Growth rate parameter; interesting values range over[0, 4].- Variables:
r (float) – Growth rate parameter.
- class physicskit.chaos.Lorenz(sigma=10.0, rho=28.0, beta=2.6666666666666665)[source]#
Bases:
DynamicalSystemThe Lorenz attractor.
- Parameters:
- Variables:
beta (sigma, rho,) – System parameters.
- property params: NDArray[float64]#
Parameter vector
(sigma, rho, beta).- Returns:
ndarray of float, shape (3,)
- class physicskit.chaos.MagneticPendulum(magnet_positions=None, friction=0.2, spring=0.2, height=0.2, strength=1.0)[source]#
Bases:
DynamicalSystemA 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.
- property params: NDArray[float64]#
Parameter vector
(friction, spring, height, strength, mx_0, my_0, ...).- Returns:
ndarray of float, shape (4 + 2*n_magnets,)
- class physicskit.chaos.QuantumBakersMap(dim, alpha=0.5)[source]#
Bases:
objectThe 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 theq < alphaandq >= alphaposition-basis blocks (mirroring the classical map stretching each piece independently), then transforms the result back to the full position representation.- Parameters:
- Variables:
- Raises:
InvalidParameterError – If alpha does not satisfy
0 < alpha < 1, or if roundingalpha * dimto 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 * dimto be an integer; for other values, the nearest integer split is used, andalphais left at the value the caller requested rather than silently adjusted to the value actually realized – pass a dim that makesalpha * dim(near-)integral for the closest match.- eigenphases()[source]#
Quasi-energies (eigenphases) of the Floquet operator.
Hand these to
physicskit.rmtto study their spacing statistics.
- evolve(psi, n_steps=1)[source]#
Propagate a state through
n_stepsmap iterations.- Parameters:
- Return type:
- Returns:
ndarray of complex, shape (n_steps + 1, dim) – The state after each iteration, including the (normalized) initial state as row 0.
- class physicskit.chaos.QuantumBilliard(billiard, resolution=150)[source]#
Bases:
objectDirichlet 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 eigenvaluesE_n = k_n^2(in units wherehbar^2 / 2m = 1) obey Weyl’s law on average (seeweyl_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
BilliardSystemalready in physicskit.chaos (no shape-specific code needed): a regular grid is laid over the shape’s bounding box, points inside the boundary (viapoints_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:
- 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:
- 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 arenan(soplt.imshow/pcolormeshleave them blank).
- 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 statisticphysicskit.rmtis built to quantify, usingwavenumbers()as its input.
- class physicskit.chaos.QuantumKickedRotor(k=1.0, dim=64, hbar=None)[source]#
Bases:
objectThe 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 classicalStandardMap(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 classicalStandardMaporbit, 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 classicalStandardMap’s k exactly.dim (
int) – Hilbert space dimension (number of angle basis states).hbar (
float|None) – Effective Planck constant; defaults to2*pi/dim, the standard choice that keeps the quantized torus’s phase-space cell count equal to dim.
- Variables:
- 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.
- eigenphases()[source]#
Quasi-energies (eigenphases) of the Floquet operator.
Hand these to
physicskit.rmtto study their spacing statistics.
- evolve(psi, n_steps=1)[source]#
Propagate a state through
n_stepskicks.- Parameters:
- Return type:
- Returns:
ndarray of complex, shape (n_steps + 1, dim) – The state after each kick, including the (normalized) initial state as row 0.
- class physicskit.chaos.RectangleBilliard(width=2.0, height=1.0)[source]#
Bases:
_RayTracingBilliardIntegrable rectangular billiard centered at the origin.
- Parameters:
- Variables:
- class physicskit.chaos.RestrictedThreeBody(mu=0.012277471)[source]#
Bases:
DynamicalSystemThe planar circular restricted three-body problem (CR3BP).
A massless third body moves under the gravity of two massive primaries (masses
1 - muandmu, 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 withphysicskit.chaos.core.integrators.rk4_integrate()(the symplecticleapfrog_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.
- jacobi_constant(state)[source]#
The Jacobi constant: the CR3BP’s conserved rotating-frame energy analog.
- Parameters:
- Return type:
- Returns:
float – The Jacobi constant
C = 2*Omega(x, y) - (vx^2 + vy^2), whereOmegais the effective (gravitational + centrifugal) potential. Conserved along any trajectory; useful for checking integrator fidelity (e.g. withphysicskit.chaos.utils.metrics.energy_drift()).
- property params: NDArray[float64]#
Parameter vector
(mu,).- Returns:
ndarray of float, shape (1,)
- class physicskit.chaos.Rossler(a=0.2, b=0.2, c=5.7)[source]#
Bases:
DynamicalSystemThe Rossler attractor.
- Parameters:
- Variables:
c (a, b,) – System parameters.
- property params: NDArray[float64]#
Parameter vector
(a, b, c).- Returns:
ndarray of float, shape (3,)
- class physicskit.chaos.SinaiBilliard(cell_size=2.0, scatterer_radius=0.5)[source]#
Bases:
_RayTracingBilliardChaotic (defocusing) Sinai billiard.
A square cell with a circular scatterer removed from its center.
- Parameters:
- Variables:
- Raises:
ValueError – If scatterer_radius is not smaller than half of cell_size.
- class physicskit.chaos.SmaleHorseshoe(contraction=0.3333333333333333, expansion=3.0)[source]#
Bases:
DiscreteMapSmale’s horseshoe map: stretch the unit square, fold it, lay it back across itself.
The square
Q = [0, 1] x [0, 1]is squeezed horizontally bylam < 1/2, stretched vertically bymu > 2, and bent into a U whose two legs lie back acrossQas the vertical stripsV0 = [0, lam] x [0, 1]andV1 = [1 - lam, 1] x [0, 1]. Only the horizontal stripsH0 = {y <= 1/mu}andH1 = {y >= 1 - 1/mu}return toQ, each by an affine branch (S. Smale, “Differentiable Dynamical Systems,” Bull. Amer. Math. Soc. 73, 747-817 (1967), Sec. I.5; J. Guckenheimer and P. Holmes, Nonlinear Oscillations, Dynamical Systems, and Bifurcations of Vector Fields, Springer (1983), Sec. 5.1):\[\begin{split}f(x, y) = \begin{cases} (\lambda x,\ \mu y) & (x, y) \in H_0 \\ (1 - \lambda x,\ \mu (1 - y)) & (x, y) \in H_1 \end{cases}\end{split}\]The second branch is orientation-reversing: that is the fold. The middle strip maps onto the bend of the U, outside
Q.The points that stay in
Qunder all forward and backward iterates form a Cantor setLambda(the product of two middle-thirds-like Cantor sets), on whichfis conjugate to the full shift on two symbols: each bi-infinite sequence ofH0/H1visits is realized by exactly one orbit. Hencef^nhas exactly2^nfixed points inQ(seeperiodic_points()) and the topological entropy isln 2.- Parameters:
- Variables:
- Raises:
ValueError – If contraction is not in
(0, 1/2)or expansion is not greater than 2.
Examples
>>> horseshoe = SmaleHorseshoe() >>> horseshoe.periodic_points(3).shape (8, 2)
- fold(points, fraction=1.0)[source]#
Image of points of the square, part-way through the stretch and fold.
- periodic_points(period)[source]#
All
2^periodfixed points off^period, in closed form.Each branch acts separately and affinely on
xandy, so the composition along an itinerarys_0 ... s_{n-1}(s_k = 1if thek-th iterate is inH1) isx -> a x + b,y -> c y + dwith fixed point(b / (1 - a), d / (1 - c)).
- survives(points, n_forward=0, n_backward=0)[source]#
Whether points stay in the square for the given numbers of iterations.
The points surviving
nforward iterations fill2^nhorizontal strips, those survivingnbackward iterations2^nvertical strips; both together,4^nsquares converging on the invariant Cantor setLambda.
- class physicskit.chaos.StandardMap(k=1.0)[source]#
Bases:
DiscreteMapThe Chirikov-Taylor standard map on the
(theta, p)cylinder.k=0is integrable; chaos onset is aroundk~1, with global chaos fork >~ 4-5.
- class physicskit.chaos.TruncatedCircleBilliard(radius=1.0, cut=0.3)[source]#
Bases:
_RayTracingBilliardA disk truncated by a straight chord.
The region
x <= radius - cutof the disk of the given radius is kept.- Parameters:
- Variables:
- Raises:
ValueError – If cut does not satisfy
0 < cut < radius.
- physicskit.chaos.kuramoto_order_parameter_lorentzian(K, gamma=1.0)[source]#
Exact steady-state order parameter for Lorentzian natural frequencies.
\[\begin{split}r(K) = \begin{cases} 0, & K \le 2\gamma, \\ \sqrt{1 - 2\gamma/K}, & K > 2\gamma. \end{cases}\end{split}\]- Parameters:
- Return type:
- Returns:
ndarray of float
Examples
>>> kuramoto_order_parameter_lorentzian([1.0, 4.0], gamma=1.0).tolist() [0.0, 0.7071067811865476]