symplectic

symplectic.py — Unified symplectic geometry toolkit for Hamiltonian mechanics

Overview

The symplectic module provides a comprehensive set of tools for studying Hamiltonian systems with an arbitrary number of degrees of freedom. It is built around the canonical symplectic structure ω = Σ dxᵢ ∧ dpᵢ and offers:

• Construction of the symplectic form and its inverse (Poisson tensor) for any n.

• Integration of Hamilton’s equations using symplectic integrators (symplectic Euler, velocity Verlet) as well as standard RK45.

• Symbolic computation of Poisson brackets.

• Detection and linear stability analysis of fixed points (equilibria).

• 1‑DOF specific utilities: action‑angle variables, phase portraits, separatrix analysis.

• 2‑DOF specific utilities: Poincaré sections, first‑return maps, monodromy matrix, Lyapunov exponents.

• Automatic inference of phase‑space variables from the Hamiltonian expression.

• Integrability classification utilities: level-spacing statistics (Berry-Tabor / BGS), Weyl’s law, Berry-Tabor smoothed density of states, KAM tori detection, and winding / rotation numbers (IntegrabilityAnalysis class).

The module is designed to be both a pedagogical tool for exploring Hamiltonian dynamics and a practical library for more advanced studies (e.g., computing monodromy matrices of periodic orbits, classifying systems as integrable or chaotic).

Mathematical background

A Hamiltonian system with n degrees of freedom is defined on a phase space of dimension 2n with canonical coordinates (x₁, p₁, …, xₙ, pₙ). The symplectic form

ω = Σ dxᵢ ∧ dpᵢ

is a closed non‑degenerate 2‑form that encodes the geometry of the space. Hamilton’s equations are derived from the relation ι_{X_H} ω = dH and read

ẋᵢ = ∂H/∂pᵢ ,  ṗᵢ = –∂H/∂xᵢ .

The Poisson bracket of two observables f, g is

{f, g} = Σ (∂f/∂xᵢ ∂g/∂pᵢ – ∂f/∂pᵢ ∂g/∂xᵢ)

and satisfies the Jacobi identity. It corresponds to the commutator in quantum mechanics via the Dirac rule {f, g} → (1/iℏ)[f̂, ĝ].

Symplectic integrators preserve the symplectic structure exactly (up to machine precision) and therefore exhibit excellent long‑term energy conservation. The module implements:

• Symplectic Euler (first‑order)

• Velocity Verlet (second‑order, equivalent to the Störmer–Verlet method)

For non‑symplectic comparisons, a standard ‘rk45’ integrator is also available.

For 1‑DOF systems, the phase space is two‑dimensional; energy surfaces are curves. The action variable

I(E) = (1/2π) ∮ p dx

is an adiabatic invariant and, for integrable systems, can be used to construct action‑angle coordinates (I, θ) in which H = H(I) and θ̇ = ω(I) = dH/dI.

For 2‑DOF systems, the phase space is four‑dimensional. A Poincaré section reduces the dynamics to a two‑dimensional map, revealing regular and chaotic motion. The monodromy matrix (linearised return map) of a periodic orbit has eigenvalues (Floquet multipliers) that characterise its stability. Lyapunov exponents quantify the rate of separation of nearby trajectories.

References

[1]

Arnold, V. I. Mathematical Methods of Classical Mechanics, Springer‑Verlag, 1989 (2nd ed.). Chapters 7–9.

[2]

Goldstein, H., Poole, C., & Safko, J. Classical Mechanics, Addison‑Wesley, 2002 (3rd ed.). Chapter 8.

[3]

Hairer, E., Lubich, C., & Wanner, G. Geometric Numerical Integration, Springer, 2006 (2nd ed.). Chapters 1, 6.

[4]

Lichtenberg, A. J., & Lieberman, M. A. Regular and Chaotic Dynamics, Springer, 1992.

[5]

Meyer, K. R., Hall, G. R., & Offin, D. Introduction to Hamiltonian Dynamical Systems and the N‑Body Problem, Springer, 2009 (2nd ed.).

[6]

Berry, M. V. & Tabor, M. “Level clustering in the regular spectrum”, Proc. R. Soc. Lond. A 356, 375–394, 1977.

[7]

Bohigas, O., Giannoni, M. J., & Schmit, C. “Characterization of chaotic quantum spectra and universality of level fluctuation laws”, Phys. Rev. Lett. 52, 1–4, 1984.

class symplectic.IntegrabilityAnalysis

Bases: object

Spectral and topological tools for classifying Hamiltonian systems as integrable, chaotic, or intermediate.

These methods are the symplectic counterparts of the geometric indicators computed in geometry.py (Poincaré sections, Lyapunov exponents, monodromy matrices). Together they provide a complete picture of the regular/chaotic nature of a system:

Spectral statistics

• Level-spacing statistics (analyze_integrability) — quantum/semiclassical fingerprint: Poisson distribution for integrable systems (Berry-Tabor conjecture, 1977), Wigner surmise for chaotic ones (BGS conjecture, 1984).

• Brody distribution (brody_distribution) — one-parameter family P(s; β) interpolating continuously between Poisson (β=0, integrable) and Wigner-GOE (β=1, chaotic). The fitted parameter β ∈ [0, 1] quantifies the degree of level repulsion and serves as a scalar chaos indicator.

Semiclassical density of states

• Weyl’s law (weyl_law) — asymptotic count of states N(E) ~ Vol{H≤E}/(2πℏ)ⁿ.

• Berry-Tabor smoothed density (berry_tabor_formula) — semiclassical density of states built from a sum over periodic orbits, each contributing a Gaussian peak weighted by its period.

Topological structure

• KAM tori detection (detect_kam_tori) — clusters periodic orbits by action proximity using Ward hierarchical clustering; each cluster corresponds to a KAM torus.

• Winding / rotation numbers (winding_number, rotation_numbers) — topological invariants that distinguish resonant (rational ratio) and quasi-periodic (irrational ratio) motion on invariant tori.

• Topological monodromy (topological_monodromy) — detects non-trivial monodromy of the action-angle fibration (Duistermaat 1980). Given a second integral of motion L, transports the action vector (I₁, I₂) around a closed loop in the energy-momentum plane (E, ℓ) encircling a critical value. The resulting integer matrix M ∈ GL(2, ℤ) is the identity for trivially fibred systems and M = [[1,1],[0,1]] for focus-focus singularities (e.g. the spherical pendulum), signalling a global obstruction to smooth action-angle coordinates.

Scarring

• Scar intensity (scar_intensity) — classical precursor of quantum scarring (Heller 1984). Measures the fraction of trajectory time spent within a prescribed neighbourhood of a periodic orbit and compares it to the ergodic baseline, yielding a scalar S ≥ 0 where S ≫ 1 indicates that the trajectory is anomalously concentrated along the orbit.

References

[BT77]

Berry, M. V. & Tabor, M., “Level clustering in the regular spectrum”, Proc. R. Soc. Lond. A 356, 375–394 (1977).

[BGS84]

Bohigas, O., Giannoni, M. J., & Schmit, C., “Characterization of chaotic quantum spectra and universality of level fluctuation laws”, Phys. Rev. Lett. 52, 1–4 (1984).

[Ar89]

Arnold, V. I., Mathematical Methods of Classical Mechanics, Springer-Verlag, 1989, Chapters 10–11.

[Br73]

Brody, T. A., “A statistical measure for the repulsion of energy levels”, Lett. Nuovo Cimento 7, 482–484 (1973).

[Du80]

Duistermaat, J. J., “On global action-angle coordinates”, Comm. Pure Appl. Math. 33, 687–706 (1980).

[Cu97]

Cushman, R. H. & Bates, L. M., Global Aspects of Classical Integrable Systems, Birkhäuser, 1997.

[He84]

Heller, E. J., “Bound-state eigenfunctions of classically chaotic Hamiltonian systems: scars of periodic orbits”, Phys. Rev. Lett. 53, 1515–1518 (1984).

static analyze_integrability(H=None, vars_phase=None, *, levels: ndarray = None, spacings: ndarray = None, traj: dict = None, ndof: int = None, second_integrals=None, lyapunov_traj: dict = None, lyapunov_dt: float = None, unfold_degree: int = 5, min_spacings: int = 30, naff_keys: tuple = None) → dict

Multi-evidence integrability classifier for Hamiltonian systems.

Design principles

1. Own the pipeline where possible. When raw energy levels are provided they are unfolded internally (polynomial fit to the staircase) before any statistic is computed, removing the system-specific smooth background that pollutes spacing tests.

2. Independent evidence channels. Five channels are evaluated and each produces its own conclusion before any aggregation:

   • Spectral — Brody β fitted by MLE plus two-sided KS tests against the Poisson and Wigner reference CDFs. These are the same data viewed through complementary lenses and are kept separate rather than averaged.

   • Algebraic — for each candidate second integral Lᵢ, the Poisson bracket {H, Lᵢ} is computed symbolically. A confirmed {H, L} = 0 is a rigorous proof of a conserved quantity and acts as a hard gate: the verdict is forced to 'Integrable' regardless of all other channels.

   • Frequency — rotation/winding numbers extracted by windowed FFT (NAFF-lite) from trajectory data, not by endpoint unwrapping. This makes the estimate robust to integration length and avoids the wrap-count artefacts of the arctan2 approach.

   • Lyapunov — the maximal Lyapunov exponent λ₁ is the cleanest dynamical chaos indicator. λ₁ > threshold is a hard chaotic gate; λ₁ ≈ 0 is a necessary (though not sufficient) condition for integrability.

   • Topological — scar intensity S is reported as a diagnostic but does not enter the verdict score because its connection to integrability is indirect (see notes).

3. Hard gates before soft scoring. The algebraic and Lyapunov channels can each override the soft score:

   • Any {H, L} = 0 with L independent of H → verdict 'Integrable' (algebraic proof, score forced to 1.0).

   • λ₁ > lyapunov_chaos_threshold → verdict 'Chaotic' (dynamical proof, score forced to 0.0).

4. Calibrated soft score for the remainder. Only the spectral channel (Brody β and KS p-values) and the frequency channel contribute to the soft score; they are combined with explicit weights reflecting their differing reliability. The ratio R is still reported for backward compatibility but does not enter the score because it is a low-power redundant summary of β.

param H:

Hamiltonian expression. Required for the algebraic channel.

type H:

sympy.Expr, optional

param vars_phase:

Phase-space variables [x₁, p₁, ...]. Required when H is provided.

type vars_phase:

list of sympy.Symbol, optional

param levels:

Raw (unsorted, unnormalised) energy eigenvalues. When provided, they are unfolded internally before computing spacing statistics. Preferred over spacings because unfolding is then correct.

type levels:

ndarray of shape (N,), optional

param spacings:

Pre-computed nearest-neighbour spacings (already unfolded to unit mean, or raw — normalisation is applied internally). Used only when levels is not given. When both are supplied levels takes precedence.

type spacings:

ndarray of shape (N,), optional

param traj:

Trajectory dict from hamiltonian_flow(). Activates the frequency channel (NAFF rotation/winding numbers).

type traj:

dict, optional

param ndof:

Number of degrees of freedom. Inferred from vars_phase when possible; required when only traj is provided.

type ndof:

int, optional

param second_integrals:

Candidate conserved quantities L₁, L₂, … Each is tested via {H, Lᵢ} = 0. Even a single confirmed integral activates the algebraic hard gate.

type second_integrals:

sympy.Expr or list of sympy.Expr, optional

param lyapunov_traj:

Separate (typically longer) trajectory used for Lyapunov estimation. Falls back to traj if not provided.

type lyapunov_traj:

dict, optional

param lyapunov_dt:

Time step for the Lyapunov QR algorithm. Inferred from lyapunov_traj['t'] when not given.

type lyapunov_dt:

float, optional

param unfold_degree:

Polynomial degree for spectral unfolding (default 5). Reduce to 3 for small samples (N < 100).

type unfold_degree:

int

param min_spacings:

Minimum number of spacings required to run spectral tests (default 30).

type min_spacings:

int

param naff_keys:

(x_key, p_key) for 1-DOF, or (x1_key, p1_key, x2_key, p2_key) for 2-DOF. If None, standard names ('x'/'p' or 'x1'/'p1'/ 'x2'/'p2') are tried automatically.

type naff_keys:

tuple of str, optional

returns:

Always present:

• verdict – 'Integrable', 'Likely integrable', 'Mixed', 'Likely chaotic', or 'Chaotic'

• verdict_source – which channel(s) determined the verdict: 'algebraic_proof', 'lyapunov_gate', or 'soft_score'

• soft_score – float in [0, 1]; 1 = integrable, 0 = chaotic. None when overridden by a hard gate.

• summary – formatted multi-line diagnostic report

• channels – dict of per-channel sub-results (see below)

• warnings – list of non-fatal diagnostic messages

channels sub-keys (present only when the channel was active):

• spectral — {'beta', 'beta_std', 'ks_poisson_p', 'ks_wigner_p', 'n_spacings', 'unfolded', 'spacings_norm'}

• algebraic — {'brackets', 'any_conserved', 'independent_integrals'} where brackets is a list of {'L': expr, 'bracket': expr, 'is_zero': bool, 'independent': bool}

• frequency — {'omega1', 'omega2', 'ratio', 'ratio_fraction', 'is_rational', 'method': 'naff'} (1-DOF: {'omega', 'method': 'naff'})

• lyapunov — {'lambda_max', 'exponents', 'is_chaotic', 'threshold'}

rtype:

dict

raises ValueError:

If neither levels nor spacings is provided AND no trajectory is given AND no symbolic data is available — i.e. if there is literally nothing to analyse.

Notes

Scar intensity is not included in the verdict. Scarring (S ≫ 1) means the trajectory is concentrated near an unstable periodic orbit, which is a property of the specific trajectory rather than the system. A scarred trajectory in a chaotic system is not evidence of integrability, and an unscarred trajectory in an integrable system is not evidence of chaos. Use scar_intensity() directly as a standalone diagnostic.

Rotation number reliability. The NAFF-lite estimator is more robust than endpoint unwrapping for trajectories of moderate length, but it still requires the trajectory to span several oscillation periods. For very short trajectories (fewer than ~10 oscillations) the frequency estimate may be unreliable; a warning is added to result['warnings'].

KAM tori detection is not included. detect_kam_tori() clusters user-supplied orbit dicts and does not detect tori from trajectory data. Its result depends entirely on what orbits the caller provides, making it unsuitable as an automatic evidence channel. Use it as a standalone tool for orbit families you have already computed.

Lyapunov convergence. The built-in lyapunov_exponents() uses a simplified QR scheme; it requires long trajectories to converge (tmax ≫ 1/λ₁). For conclusive Lyapunov evidence integrate for at least 10³ natural periods.

References

[BT77]

Berry & Tabor, Proc. R. Soc. Lond. A 356 (1977).

[BGS84]

Bohigas, Giannoni & Schmit, PRL 52 (1984).

[Br73]

Brody, Lett. Nuovo Cimento 7 (1973).

[Las90]

Laskar, Icarus 88 (1990) — NAFF algorithm.

[OC92]

Ott, Chaos in Dynamical Systems, Cambridge (1992).

Examples

Integrable — symbolic proof (no spectrum needed):

>>> x1, p1, x2, p2 = symbols('x1 p1 x2 p2', real=True)
>>> H = (p1**2 + x1**2 + p2**2 + x2**2) / 2
>>> L = (p2**2 + x2**2) / 2
>>> r = IntegrabilityAnalysis.analyze_integrability(
...         H=H, vars_phase=[x1,p1,x2,p2], second_integrals=L)
>>> r['verdict']
'Integrable'
>>> r['verdict_source']
'algebraic_proof'

From a raw spectrum:

>>> levels = np.sort(rng.exponential(1.0, 500).cumsum())
>>> r = IntegrabilityAnalysis.analyze_integrability(levels=levels)
>>> r['channels']['spectral']['beta']   # close to 0
< 0.2

From a trajectory (frequency channel):

>>> traj = hamiltonian_flow(H, (1,0,0,1), (0, 60*np.pi),
...                         vars_phase=[x1,p1,x2,p2], n_steps=8000)
>>> r = IntegrabilityAnalysis.analyze_integrability(
...         traj=traj, ndof=2,
...         H=H, vars_phase=[x1,p1,x2,p2], second_integrals=L)
>>> r['verdict']
'Integrable'

static berry_tabor_formula(orbits, energy: float, window: float = 1.0) → float

Berry-Tabor smoothed density of states from a collection of periodic orbits.

Each orbit contributes a Gaussian peak centred at its energy:

ρ(E) = Σ_γ T_γ · exp(-(E_γ - E)² / 2σ²) / (σ√(2π) · 2π)

where T_γ is the period and σ = window is the energy smoothing width.

Parameters:

• orbits (iterable) – Each element must expose .energy and .period attributes or be a dict with 'energy' and 'period' keys (as returned, e.g., by the monodromy / periodic-orbit routines in this module).

• energy (float) – Evaluation energy E.

• window (float) – Gaussian smoothing width σ in energy units (default 1.0).

Returns:

Smoothed density of states ρ(E).

Return type:

float

Examples

>>> # Using plain dicts (compatible with symplectic.py orbit dicts)
>>> orbits = [{'energy': 1.0, 'period': 6.28},
...           {'energy': 2.0, 'period': 6.28}]
>>> IntegrabilityAnalysis.berry_tabor_formula(orbits, energy=1.5)
...

Notes

For integrable systems the smoothed density converges to the Weyl term plus Berry-Tabor oscillatory corrections. For chaotic systems use the Gutzwiller trace formula instead (implemented in geometry.py).

static brody_distribution(spacings: ndarray, fit: bool = True) → dict

Brody’s intermediate level-spacing distribution.

Brody (1973) introduced a one-parameter family that interpolates continuously between the Poisson distribution (integrable, β=0) and the Wigner surmise (chaotic/GOE, β=1):

P(s; β) = (β+1) b s^β exp(-b s^{β+1})

where b = [Γ((β+2)/(β+1))]^{β+1} ensures unit mean spacing.

Parameters:

• spacings (ndarray of shape (N,)) – Raw nearest-neighbour energy-level spacings (need not be pre-normalised; normalisation is performed internally).

• fit (bool) – If True (default), estimate β from the data by maximum likelihood. If False, return only the theoretical curves for β=0 and β=1 without fitting.

Returns:

• beta – fitted Brody parameter (0 ≤ β ≤ 1),

  or None if fit=False

• beta_std – standard error of the fit (bootstrap),

  or None if fit=False

• classification – ‘Integrable (β≈0)’, ‘Chaotic (β≈1)’,

  or ‘Intermediate (β=<value>)’

• pdf – callable P(s; β_fit) for plotting

• nll – negative log-likelihood at optimum

Return type:

dict with keys

Examples

>>> rng = np.random.default_rng(42)
>>> wigner_spacings = rng.rayleigh(scale=np.sqrt(4/np.pi), size=500)
>>> result = IntegrabilityAnalysis.brody_distribution(wigner_spacings)
>>> round(result['beta'], 1)
1.0

Notes

The fit uses scipy.optimize.minimize_scalar on the negative log-likelihood, which is unimodal in β ∈ [0,1] for typical spectra. Bootstrap standard error uses 200 resamples.

References

[Br73]

Brody, T. A., “A statistical measure for the repulsion of energy levels”, Lett. Nuovo Cimento 7, 482–484 (1973).

static detect_kam_tori(orbits, tolerance: float = 0.1) → dict

Cluster periodic orbits into KAM tori by action proximity.

In an integrable (or near-integrable) system, periodic orbits of nearby action lie on the same KAM torus. This function uses Ward hierarchical clustering on the action values to group orbits accordingly.

Parameters:

• orbits (iterable) – Each element must expose .action, .energy, .period and .stability / 'stability' (or 'stability_1' for 2-DOF objects) or be a dict with the same keys.

• tolerance (float) – Maximum within-cluster action distance for the Ward linkage cut (default 0.1).

Returns:

• n_tori – number of detected tori

• tori – list of dicts, one per torus, each containing:

  • id : cluster id (int)

  • n_orbits : number of member orbits

  • action : mean action of member orbits

  • energy : mean energy of member orbits

  • period : mean period of member orbits

  • stable : True if mean stability exponent < 0

Return type:

dict with keys

Examples

>>> result = IntegrabilityAnalysis.detect_kam_tori(orbits, tolerance=0.05)
>>> print(result['n_tori'], result['tori'][0]['action'])

Notes

At least two orbits are required for clustering; a single orbit is returned as a degenerate torus.

static rotation_numbers(traj: dict, x_key: str = 'x1', p_key: str = 'p1', y_key: str = 'x2', q_key: str = 'p2') → tuple

Rotation numbers (ω_1, ω_2) of a 2-DOF trajectory.

Each rotation number is the average angular velocity in the corresponding (xᵢ, pᵢ) phase-space plane:

ωᵢ = (θᵢ(T) - θᵢ(0)) / (2π · T), θᵢ = arctan2(pᵢ, xᵢ).

Their ratio ω_1/ω_2 is the frequency ratio; rational values indicate resonant (periodic) orbits, irrational values indicate quasi-periodic motion on a KAM torus.

Parameters:

• traj (dict) – Trajectory dict as returned by hamiltonian_flow.

• x_key (str) – Keys for the first position and momentum (default 'x1', 'p1').

• p_key (str) – Keys for the first position and momentum (default 'x1', 'p1').

• y_key (str) – Keys for the second position and momentum (default 'x2', 'p2').

• q_key (str) – Keys for the second position and momentum (default 'x2', 'p2').

Returns:

Rotation numbers (ω_1, ω_2).

Return type:

(float, float)

Examples

>>> x1,p1,x2,p2 = symbols('x1 p1 x2 p2', real=True)
>>> H = (p1**2 + p2**2 + x1**2 + x2**2) / 2   # isotropic oscillator
>>> traj = hamiltonian_flow(H, (1,0,0,1), (0, 20*np.pi),
...                         vars_phase=[x1,p1,x2,p2], n_steps=5000)
>>> IntegrabilityAnalysis.rotation_numbers(traj)
(0.5, 0.5)   # both oscillators at frequency 1/2π

Notes

For short trajectories the rotation number estimate may be noisy; use longer integration times for robust results.

static scar_intensity(traj: dict, vars_phase, orbit_points: ndarray, radius: float = None) → dict

Measure the scar intensity of a trajectory relative to a periodic orbit.

Scarring (Heller 1984) is the tendency of trajectories near an unstable periodic orbit to spend more time there than ergodicity would predict. The classical signature is a dwell-time excess: the fraction of time steps within distance radius of the orbit is larger than the fraction of phase-space volume that ball occupies.

The scar intensity S is defined as:

S = f_orbit / f_expected

where:

f_orbit = fraction of trajectory points within distance r of

any orbit point

f_expected = π r² / A_traj (area of the ball relative to the

bounding box of the trajectory — a crude ergodic baseline)

S ≈ 1 → ergodic (no scarring) S >> 1 → trajectory is scarred along the orbit

Parameters:

• traj (dict) – Trajectory dict as returned by hamiltonian_flow. For 2-DOF systems the projection onto (x1, p1) is used.

• vars_phase (list of sympy.Symbol) – [x, p] or [x1, p1, x2, p2].

• orbit_points ((K, 2) ndarray) – Points (x, p) sampled along the periodic orbit.

• radius (float, optional) – Proximity radius r for dwell-time counting. Defaults to 1/10 of the trajectory’s diameter in phase space.

Returns:

• scar_intensity – S = f_orbit / f_expected (float)

• f_orbit – fraction of trajectory within r of orbit

• f_expected – ergodic baseline fraction

• n_close – number of trajectory points within r

• radius – radius used

Return type:

dict with keys

Examples

>>> x, p = symbols('x p', real=True)
>>> H = (p**2 + x**2) / 2
>>> traj = hamiltonian_flow(H, (1, 0), (0, 40*np.pi),
...                         vars_phase=[x, p], n_steps=4000)
>>> theta = np.linspace(0, 2*np.pi, 100, endpoint=False)
>>> orbit = np.column_stack([np.cos(theta), np.sin(theta)])
>>> result = IntegrabilityAnalysis.scar_intensity(traj, [x, p], orbit)
>>> result['scar_intensity']   # >> 1 since traj IS the orbit
> 5.0

References

[He84]

Heller, E. J., “Bound-state eigenfunctions of classically chaotic Hamiltonian systems: scars of periodic orbits”, Phys. Rev. Lett. 53, 1515–1518 (1984).

static topological_monodromy(H, L, vars_phase, critical_value: tuple, loop_radius: float = 0.5, n_loop_points: int = 48) → dict

Detect non-trivial topological monodromy (Duistermaat 1980).

For a 2-DOF integrable system with Hamiltonians H and a second integral L, the energy-momentum map F: M → ℝ² sends each phase-space point to (H, L). The fibres of F are Liouville tori parametrised by action variables (I₁, I₂).

Around a critical value (E*, ℓ*) of F (a focus-focus singularity), the action lattice undergoes a non-trivial linear transport:

(I₁, I₂) → M · (I₁, I₂), M ∈ GL(2, ℤ)

For a focus-focus point M = [[1, 1], [0, 1]] (upper-triangular, non-identity), signalling a global obstruction to action-angle variables (Duistermaat 1980; Cushman & Bates 1997).

The loop is traced in the (E, ℓ) plane as:

E(θ) = E* + r·cos(θ), ℓ(θ) = ℓ* + r·sin(θ), θ ∈ [0, 2π).

At each (E(θ), ℓ(θ)) the two actions are computed from the 1-DOF sub-problems obtained by energy-momentum reduction, and the transport matrix is read off at the end of the loop.

Parameters:

• H (sympy.Expr) – 2-DOF Hamiltonian H(x1, p1, x2, p2).

• L (sympy.Expr) – Second conserved quantity (must Poisson-commute with H).

• vars_phase (list of sympy.Symbol) – [x1, p1, x2, p2].

• critical_value ((float, float)) – (E*, ℓ*) — the critical value of the energy-momentum map around which the loop is traced.

• loop_radius (float) – Radius r of the loop in (E, ℓ) space (default 0.5).

• n_loop_points (int) – Number of sample points along the loop (default 48).

Returns:

• monodromy_matrix – 2×2 integer ndarray M (rounded from float)

• monodromy_float – 2×2 float ndarray before rounding

• is_trivial – True if M == identity

• actions_along_loop – (n_loop_points, 2) array of (I₁(θ), I₂(θ))

• angles – loop parameter θ ∈ [0, 2π)

• loop_EL – (n_loop_points, 2) array of (E(θ), ℓ(θ))

Return type:

dict with keys

Raises:

ValueError – If vars_phase is not 2-DOF.

Notes

The monodromy matrix is estimated as:

M_float = I_end · pinv(I_start_matrix)

where I_start_matrix is a 2×2 matrix built from the action vectors at θ=0 and θ=π/2 (two independent directions), giving the full linear map. Rounding to integers reflects the lattice structure of the torus.

A well-resolved loop (n_loop_points ≥ 36, loop_radius small enough to contain only one critical value) is needed for accurate results.

References

[Du80]

Duistermaat, J. J., “On global action-angle coordinates”, Comm. Pure Appl. Math. 33, 687–706 (1980).

[Cu97]

Cushman, R. H. & Bates, L. M., Global Aspects of Classical Integrable Systems, Birkhäuser, 1997.

static unfold_spectrum(levels: ndarray, poly_degree: int = 5) → ndarray

Unfold a raw energy spectrum to unit mean level density.

Unfolding maps eigenvalues E_i to dimensionless levels ε_i so that the mean nearest-neighbour spacing is exactly 1, removing the system-specific smooth part of the density of states (the Weyl term) and leaving only the universal fluctuation statistics.

The smooth cumulative level count N̄(E) is estimated by fitting a polynomial of degree poly_degree to the staircase function N(E) = #{levels ≤ E}. The unfolded levels are then ε_i = N̄(E_i).

Parameters:

• levels (ndarray of shape (N,)) – Raw sorted (or unsorted) energy eigenvalues.

• poly_degree (int) – Degree of the polynomial used to fit the smooth spectral staircase N̄(E). Higher degrees track non-uniform densities better but may overfit for small samples (default 5).

Returns:

Unfolded levels ε_i with unit mean spacing.

Return type:

ndarray of shape (N,)

Notes

For small samples (N < 50) reduce poly_degree to 3 to avoid overfitting. For systems whose density of states is known analytically (e.g. harmonic oscillators via Weyl’s law) prefer passing the analytic N̄ as an interpolation table rather than using this polynomial fit.

Examples

>>> levels = np.sort(np.random.default_rng(0).normal(0, 1, 200))
>>> unfolded = IntegrabilityAnalysis.unfold_spectrum(levels)
>>> spacings = np.diff(unfolded)
>>> np.isclose(spacings.mean(), 1.0, atol=0.05)
True

static weyl_law(energy: float, ndof: int, hbar: float = 1.0) → float

Weyl’s asymptotic law: number of quantum states below energy E.

The simplified form assumes that the phase-space volume scales as E^n (valid, e.g., for harmonic-oscillator-like Hamiltonians):

N(E) ≈ (1 / (2πℏ))^n · E^n

For a general Hamiltonian use phase_space_volume from this module to compute Vol{H ≤ E} and then divide by (2πℏ)^n.

Parameters:

• energy (float) – Energy threshold E.

• ndof (int) – Number of degrees of freedom n.

• hbar (float) – Reduced Planck constant (default 1.0).

Returns:

Approximate number of states N(E).

Return type:

float

Examples

>>> IntegrabilityAnalysis.weyl_law(3.0, ndof=1, hbar=1.0)
0.477...

static winding_number(traj: dict, x_key: str = 'x', y_key: str = 'y') → float

Winding number of a 2-DOF trajectory around the origin in configuration space.

The winding number counts how many times the spatial projection (x(t), y(t)) winds around the origin:

W = (θ(T) - θ(0)) / 2π, θ = arctan2(y, x).

Parameters:

• traj (dict) – Trajectory dict as returned by hamiltonian_flow.

• x_key (str) – Keys for the two configuration-space coordinates (default 'x', 'y'). Adjust if your variable names differ (e.g. 'x1', 'x2').

• y_key (str) – Keys for the two configuration-space coordinates (default 'x', 'y'). Adjust if your variable names differ (e.g. 'x1', 'x2').

Returns:

Winding number (positive = counter-clockwise).

Return type:

float

class symplectic.SymplecticForm(n=None, vars_phase=None, omega_matrix=None)

Bases: object

Represents the symplectic structure on a 2n-dimensional phase space.

The symplectic form ω is a closed, non-degenerate 2-form that defines the geometric structure of Hamiltonian mechanics. In canonical coordinates, it takes the standard form:

ω = Σᵢ dxᵢ ∧ dpᵢ

The matrix representation J (also called the symplectic matrix) satisfies: - Jᵀ = -J (antisymmetric) - J² = -I (for canonical form) - det(J) = 1

This matrix is used to convert gradients of the Hamiltonian into the Hamiltonian vector field: ż = J · ∇H

Parameters:

• n (int, optional) – Number of degrees of freedom. Used to generate generic variables if vars_phase is not provided.

• vars_phase (list of sympy.Symbol, optional) – Ordered list of canonical variables [x₁, p₁, x₂, p₂, …]. If provided, n is inferred from the length of this list.

• omega_matrix (sympy.Matrix, optional) – Custom 2n×2n symplectic matrix. If None, the canonical form is used: J = [[0, -I], [I, 0]] where I is the n×n identity.

n

Number of degrees of freedom.

Type:

int

vars_phase

Phase space variables.

Type:

list

omega_matrix

The symplectic matrix J.

Type:

sympy.Matrix

omega_inv

The inverse J⁻¹ (Poisson tensor), used for computing Poisson brackets.

Type:

sympy.Matrix

Raises:

ValueError – If neither n nor vars_phase is provided, if omega_matrix has incorrect dimensions, or if omega_matrix is not antisymmetric.

Examples

>>> # 1-DOF canonical form
>>> sf = SymplecticForm(n=1)
>>> sf.omega_matrix
Matrix([[0, -1], [1, 0]])

>>> # 2-DOF with explicit variables
>>> x1, p1, x2, p2 = symbols('x1 p1 x2 p2', real=True)
>>> sf = SymplecticForm(vars_phase=[x1, p1, x2, p2])
>>> sf.n
2

See also

poisson_bracket

Uses omega_inv to compute {f, g}

hamiltonian_flow

Uses the symplectic structure for integration

eval(point)

Evaluate the (constant) symplectic matrix at a point. For canonical form, it’s independent of coordinates.

Parameters:

point (array_like) – Phase space point (length 2n).

Returns:

2n×2n matrix ωᵢⱼ.

Return type:

ndarray

symplectic.SymplecticForm1D(vars_phase=None)

symplectic.SymplecticForm2D(vars_phase=None)

symplectic.action_angle_transform(H, x_range, p_range, vars_phase=None, n_contours=10)

Compute action-angle transformation for integrable system (1‑DOF).

For integrable 1D systems, finds action I(E) and angle θ such that:

H = H(I) θ̇ = ω(I) = dH/dI

Parameters:

• H (sympy expression) – Hamiltonian H(x, p).

• x_range (tuple) – Phase space domain.

• p_range (tuple) – Phase space domain.

• vars_phase (list, optional) – Phase space variables [x, p].

• n_contours (int) – Number of energy levels to compute.

Returns:

Contains ‘energies’, ‘actions’, ‘frequencies’.

Return type:

dict

symplectic.action_integral(H, E, vars_phase=None, method='numerical', x_bounds=None)

Compute action integral I(E) for a 1‑DOF system.

I(E) = (1/2π) ∮ p dx

Parameters:

• H (sympy expression) – Hamiltonian H(x, p).

• E (float or sympy expression) – Energy level.

• vars_phase (list, optional) – Phase space variables [x, p].

• method ({'numerical', 'symbolic'}) – Integration method.

• x_bounds (tuple, optional) – Integration limits (x_min, x_max) for the orbit.

Returns:

Action I(E).

Return type:

float or sympy expression

symplectic.evolve_phase_space_region(H, initial_region, t_eval, vars_phase=None, integrator='verlet', n_steps=10000, plot=True, ax=None, **plot_kwargs)

Evolve a closed region in phase space under the Hamiltonian flow and optionally visualise its deformation and area conservation.

The function integrates each vertex of the polygon defined by initial_region forward to the times specified in t_eval, computes the enclosed area at each time (using the shoelace formula), and plots the evolved regions if requested.

This provides a direct visual demonstration of Liouville’s theorem: the symplectic area is preserved by the Hamiltonian flow (up to numerical errors, especially when using a symplectic integrator).

Parameters:

• H (sympy.Expr) – Hamiltonian for a 1‑degree‑of‑freedom system H(x, p).

• initial_region ((N,2) array_like) – Ordered points defining a closed polygon (first and last point should be identical). The region is typically a simple shape such as a rectangle or an ellipse.

• t_eval (float or list of float) – Times at which to record the evolved region. If a single float is given, it is treated as the final time and the region is evaluated at t=0 and that time. Use a list to obtain intermediate snapshots.

• vars_phase (list of sympy.Symbol, optional) – Phase space variables [x, p]. If not provided, they are inferred from H.

• integrator ({'symplectic', 'verlet', 'rk45'}, optional) – Numerical integrator passed to hamiltonian_flow (default ‘verlet’). Symplectic integrators are recommended for accurate long‑term area preservation.

• n_steps (int, optional) – Number of integration steps used by hamiltonian_flow over the whole time interval [0, max(t_eval)]. A high value ensures accurate point positions at the requested times (default 10000).

• plot (bool, optional) – If True, generate a plot showing the region at each time in t_eval (default True).

• ax (matplotlib.axes.Axes, optional) – Axes on which to draw the plot. If None and plot=True, a new figure and axes are created.

• **plot_kwargs – Additional keyword arguments passed to the plotting call (e.g., colors, alpha, linewidth). Useful for customising the appearance of the evolved regions.

Returns:

A dictionary with the following keys: - ‘times’ : ndarray

The evaluation times (sorted, unique).

• ’region_at_t’list of (N,2) ndarrays

  Evolved polygon for each time, in the same vertex order as input.

• ’areas’ndarray

  Area enclosed by the polygon at each time (computed with the shoelace formula).

• ’plot_handles’list, optional

  If plot=True, a list of handles to the plotted lines/fills for each time (useful for further customisation or legends).

Return type:

dict

Raises:

ValueError – If the system is not 1‑DOF, if the initial region has fewer than 3 points, or if the polygon is not closed (first and last point differ beyond a small tolerance).

Examples

>>> x, p = symbols('x p', real=True)
>>> H = p**2/2 + x**2/2          # harmonic oscillator
>>> region = rectangle_region(center=(0, 1), width=0.5, height=0.3)
>>> result = evolve_phase_space_region(H, region, t_eval=[0, 2*np.pi, 4*np.pi],
...                                    integrator='verlet')
>>> result['areas']
array([0.15, 0.15, 0.15])        # area conserved (up to numerical errors)

See also

rectangle_region

Convenience function to generate a rectangular region.

phase_portrait

Plot energy contours and vector fields.

hamiltonian_flow

Underlying numerical integration routine.

symplectic.find_fixed_points(H, vars_phase=None, domain=None, tol=1e-06, numerical=False)

Find fixed points (equilibria) of Hamiltonian system: ∂H/∂z_i = 0 for all i.

Parameters:

• H (sympy expression) – Hamiltonian.

• vars_phase (list of symbols, optional) – Phase space variables.

• domain (list of tuples, optional) – Bounds for each variable [(min, max), …] for numerical search.

• tol (float) – Tolerance for numerical solutions.

• numerical (bool) – If True, use numerical root‑finding even if symbolic solve is possible.

Returns:

Fixed points (z₁, z₂, …, z₂ₙ).

Return type:

list of tuples

symplectic.first_return_map(section_points, plot_variables=('x1', 'p1'))

Construct the first return map (Poincaré map) from a Poincaré section.

For a 2-DOF Hamiltonian system, a Poincaré section reduces the 4D flow to a 2D map by recording successive intersections with a chosen surface. The first return map plots each intersection point (in the chosen coordinates) against the next one. Regular (KAM) tori appear as closed curves, while chaotic motion fills area in the return map.

This function extracts successive pairs from the list of section points (returned by poincare_section) and organizes them into two arrays: current (the i-th point) and next (the (i+1)-th point). The resulting dictionary can be directly used with plt.plot to visualise the map.

Parameters:

• section_points (list of dict) – List of section crossing points as returned by poincare_section (the ‘section_points’ field). Each dict has keys like ‘x1’, ‘p1’, ‘x2’, ‘p2’ containing the coordinates at the crossing.

• plot_variables (tuple of str, optional) – Two variable names to use as the coordinates for the return map, e.g. (‘x1’, ‘p1’) or (‘x2’, ‘p2’). Defaults to (‘x1’, ‘p1’).

Returns:

A dictionary with the following keys:

• ’current’ : ndarray, shape (N-1, 2) Coordinates of the first N-1 section points in the chosen variables.

• ’next’ : ndarray, shape (N-1, 2) Coordinates of the subsequent (i+1)-th points.

• ’variables’ : tuple of str The variable names used (same as plot_variables).

Return type:

dict

Raises:

ValueError – If fewer than 2 section points are provided.

Examples

>>> # Compute a Poincaré section
>>> ps = poincare_section(H, Sigma_def, z0, tmax=100)
>>> # Extract first return map in (x1, p1)
>>> ret_map = first_return_map(ps['section_points'], plot_variables=('x1', 'p1'))
>>> # Plot the map
>>> plt.plot(ret_map['current'][:, 0], ret_map['current'][:, 1], 'o', markersize=2)
>>> plt.xlabel(ret_map['variables'][0])
>>> plt.ylabel(ret_map['variables'][1])
>>> plt.title('First return map')
>>> plt.show()

See also

poincare_section

Compute the Poincaré section that produces the points.

visualize_poincare_section

High-level visualisation of multiple sections.

symplectic.frequency(H, I_val, method='derivative')

Compute frequency ω(I) = dH/dI from action variable.

Parameters:

• H (sympy expression) – Hamiltonian as function of action H(I).

• I_val (float or sympy expression) – Action variable value.

• method (str) – Computation method: ‘derivative’, ‘period’.

Returns:

Frequency ω(I).

Return type:

float or sympy expression

symplectic.hamiltonian_flow(H, z0, tspan, vars_phase=None, integrator='symplectic', n_steps=1000)

Numerically integrate Hamilton’s equations of motion.

Solves the system:

ẋᵢ = ∂H/∂pᵢ ṗᵢ = -∂H/∂xᵢ

for i = 1, …, n degrees of freedom.

Integrator options and their properties: - ‘symplectic’ (Symplectic Euler): First-order, exactly preserves symplectic

structure, good long-term energy behavior, computationally efficient.

• ‘verlet’ (Velocity Verlet): Second-order, time-reversible, excellent for oscillatory systems, preserves symplectic structure.

• ‘rk45’ (Runge-Kutta 4/5): Fourth-order adaptive, not symplectic, better short-term accuracy but may drift in energy over long simulations.

Parameters:

• H (sympy.Expr) – Hamiltonian function H(x₁, p₁, …, xₙ, pₙ).

• z0 (array_like) – Initial conditions as [x₁₀, p₁₀, x₂₀, p₂₀, …], length 2n.

• tspan (tuple of float) – Time integration interval (t_start, t_end).

• vars_phase (list of sympy.Symbol, optional) – Phase space variables in canonical order. If None, inferred from H.

• integrator ({'symplectic', 'verlet', 'rk45'}) – Numerical integration method (default: ‘symplectic’).

• n_steps (int) – Number of discrete time steps (default: 1000).

Returns:

Dictionary containing: - ‘t’ : ndarray, time points - ‘<var_name>’ : ndarray, trajectory for each phase space variable - ‘energy’ : ndarray, H evaluated along the trajectory

Return type:

dict

Raises:

ValueError – If integrator is not recognized, if variables cannot be inferred, or if z0 has incorrect dimension.

Notes

Performance: derivatives are lambdified once into a single vectorised function F(z) → dz that operates on a numpy state vector. The symplectic and verlet loops are written with direct array indexing (no Python-level list allocation per step), giving 10–50× speedup over the original scalar-loop implementation for large n_steps.

For long-time simulations prefer ‘symplectic’ or ‘verlet’ to preserve the geometric structure of phase space. Use ‘rk45’ for short-time high-accuracy requirements.

symplectic.hamiltonian_flow_4d(H, z0, tspan, integrator='symplectic', n_steps=1000)

Backward compatibility wrapper for 4D flow.

symplectic.linearize_at_fixed_point(H, z0, vars_phase=None)

Compute linear stability matrix (Jacobian of the Hamiltonian vector field) at a fixed point.

For Hamiltonian systems, the linearized equations are dδz/dt = J · Hessian(H) · δz, where J is the symplectic matrix. This function returns the matrix J·Hessian(H) and its eigenvalues.

Parameters:

• H (sympy expression) – Hamiltonian.

• z0 (tuple or list) – Fixed point coordinates (length 2n).

• vars_phase (list of symbols, optional) – Phase space variables.

Returns:

Contains ‘jacobian’ (numpy array), ‘eigenvalues’, ‘eigenvectors’, and a simple classification based on eigenvalues: - ‘elliptic’ if all eigenvalues are purely imaginary - ‘hyperbolic’ if all eigenvalues are real (or come in opposite pairs) - ‘mixed’ otherwise.

Return type:

dict

symplectic.lyapunov_exponents(trajectory, dt, H=None, vars_phase=None, n_vectors=None, renorm_interval=10, epsilon=1e-06)

Estimate Lyapunov exponents from a trajectory using QR algorithm.

If H and vars_phase are provided, the Jacobian of the Hamiltonian vector field is computed by finite differences, and the linearized flow is approximated as I + J_f * dt. This yields accurate exponents for integrable systems (λ≈0). Without H, a crude (and often wrong) fallback is used – a warning is issued.

Parameters:

• trajectory (dict) – Trajectory dict as returned by hamiltonian_flow.

• dt (float) – Time step between trajectory points.

• H (sympy.Expr, optional) – Hamiltonian expression (required for accurate Jacobian).

• vars_phase (list of sympy.Symbol, optional) – Phase space variables (required when H is given).

• n_vectors (int, optional) – Number of tangent vectors (default = number of DOF).

• renorm_interval (int) – Number of steps between QR renormalisations (default 10).

• epsilon (float) – Finite‑difference step for Jacobian (default 1e-6).

Returns:

Lyapunov exponents sorted in descending order.

Return type:

ndarray

Raises:

ValueError – If ndof cannot be determined.

symplectic.monodromy_matrix(H, periodic_orbit, vars_phase=None, method='finite_difference')

Compute the monodromy matrix for a periodic orbit in a 2-DOF system.

The monodromy matrix M is the linearized Poincaré return map around a periodic orbit. It describes how infinitesimal perturbations evolve over one period T:

δz(T) = M · δz(0)

The eigenvalues of M (Floquet multipliers) determine orbital stability:

• |λ| = 1 for all λ: neutrally stable (typical for Hamiltonian systems)

• Any |λ| > 1: unstable orbit

• All |λ| < 1: stable (not possible in conservative Hamiltonian systems)

For Hamiltonian systems, eigenvalues come in reciprocal pairs (λ, 1/λ) and/or complex conjugate pairs on the unit circle.

Parameters:

• H (sympy.Expr) – Hamiltonian for a 2-degree-of-freedom system.

• periodic_orbit (dict) – A trajectory dictionary (from hamiltonian_flow) representing one complete period, with keys ‘t’ and variable names.

• vars_phase (list of sympy.Symbol, optional) – Variables [x1, p1, x2, p2]. If None, uses default canonical names.

• method ({'finite_difference'}) – Computation method (default: ‘finite_difference’). Currently only finite difference perturbations are implemented.

Returns:

Dictionary containing: - ‘M’ : ndarray, 4×4 monodromy matrix - ‘eigenvalues’ : ndarray, Floquet multipliers - ‘floquet_multipliers’ : ndarray, same as eigenvalues - ‘stable’ : bool, True if all |λ| ≤ 1 (within tolerance)

Return type:

dict

Raises:

• NotImplementedError – If method is not ‘finite_difference’.

• ValueError – If the system is not 2-DOF.

Examples

>>> x1, p1, x2, p2 = symbols('x1 p1 x2 p2', real=True)
>>> H = (p1**2 + p2**2 + x1**2 + x2**2)/2
>>> orbit = hamiltonian_flow(H, [1, 0, 0, 0], (0, 2*pi),
...                          vars_phase=[x1, p1, x2, p2])
>>> mono = monodromy_matrix(H, orbit, vars_phase=[x1, p1, x2, p2])
>>> np.allclose(np.abs(mono['eigenvalues']), 1.0, atol=1e-3)
True

Notes

The finite difference method perturbs each coordinate by ε = 1e-6 and integrates forward one period. For higher accuracy, consider implementing variational equations directly.

See also

lyapunov_exponents

Long-term stability characterization

linearize_at_fixed_point

Stability analysis for equilibria

symplectic.phase_portrait(H, x_range, p_range, vars_phase=None, resolution=50, levels=20)

Generate a 2D phase portrait for a 1-DOF Hamiltonian system.

Creates a visualization showing: 1. Energy contour lines (level sets of H) 2. Hamiltonian vector field arrows indicating flow direction

For 1-DOF systems, trajectories lie on energy contours, making this an effective tool for identifying:

• Fixed points (where vector field vanishes)

• Periodic orbits (closed contours)

• Separatrices (boundaries between different types of motion)

• Allowed vs. forbidden regions

Parameters:

• H (sympy.Expr) – Hamiltonian H(x, p) for a 1-degree-of-freedom system.

• x_range (tuple of float) – Position axis limits (x_min, x_max).

• p_range (tuple of float) – Momentum axis limits (p_min, p_max).

• vars_phase (list of sympy.Symbol, optional) – Variables [x, p]. If None, inferred from H.

• resolution (int) – Grid resolution for contour and vector field computation (default: 50). Higher values give smoother plots but increase computation time.

• levels (int) – Number of energy contour levels to display (default: 20).

Returns:

Displays a matplotlib figure with the phase portrait.

Return type:

None

Raises:

ValueError – If the system is not 1-DOF or variables cannot be inferred.

Examples

>>> x, p = symbols('x p', real=True)
>>> H = p**2/2 + x**2/2  # Harmonic oscillator
>>> phase_portrait(H, x_range=(-2, 2), p_range=(-2, 2))

>>> # Pendulum
>>> H = p**2/2 - cos(x)
>>> phase_portrait(H, x_range=(-2*pi, 2*pi), p_range=(-3, 3))

Notes

The vector field is subsampled (every resolution//20 points) to avoid visual clutter. Energy contours are labeled with their H values.

See also

visualize_phase_space_structure

More comprehensive visualization with fixed points and separatrices

separatrix_analysis

Detailed analysis near saddle points

symplectic.poincare_section(H, Sigma_def, z0, tmax, vars_phase=None, n_returns=1000, integrator='symplectic')

Compute a Poincaré section for a 2-DOF Hamiltonian system.

A Poincaré section reduces the 4D continuous flow to a 2D discrete map by recording intersections of trajectories with a specified surface Σ in phase space. This reveals the underlying structure of the dynamics:

• Regular motion appears as closed curves (KAM tori)

• Chaotic motion appears as scattered points

• Periodic orbits appear as fixed points or finite sets

The section surface is defined by Σ = {(x₁, p₁, x₂, p₂) : variable = value} with an optional crossing direction.

Parameters:

• H (sympy.Expr) – Hamiltonian for a 2-degree-of-freedom system.

• Sigma_def (dict) –

  Section surface definition with keys: - ‘variable’ : str, name of the section variable (e.g., ‘x2’, ‘p1’) - ‘value’ : float, the constant value defining the surface - ‘direction’ : str, optional, ‘positive’ or ‘negative’ crossing

  (default: ‘positive’)

• z0 (array_like) – Initial condition [x₁₀, p₁₀, x₂₀, p₂₀].

• tmax (float) – Maximum integration time.

• vars_phase (list of sympy.Symbol, optional) – Variables [x1, p1, x2, p2]. If None, uses default canonical names.

• n_returns (int) – Maximum number of section crossings to collect (default: 1000).

• integrator ({'symplectic', 'verlet', 'rk45'}) – Integration method (default: ‘symplectic’).

Returns:

Dictionary containing: - ‘t_crossings’ : ndarray, times at which the section was crossed - ‘section_points’ : list of dict, phase space coordinates at each

crossing (interpolated to the exact crossing time)

Return type:

dict

Raises:

ValueError – If the system is not 2-DOF, if the section variable is not found, or if insufficient crossings are detected.

Notes

Crossing detection is fully vectorised with NumPy (no Python-level loop over time steps), giving a significant speedup for large n_steps. Crossing times are linearly interpolated between integration steps for accuracy. Use high n_steps in the underlying integration for precise section points. For chaotic systems, increase tmax and n_returns.

See also

visualize_poincare_section

Plot multiple sections together

first_return_map

Analyze the discrete map from section points

symplectic.poisson_bracket(f, g, vars_phase=None)

Compute Poisson bracket {f, g} = (∂f/∂z) · ω⁻¹ · (∂g/∂z).

Parameters:

• f (sympy expressions) – Functions on phase space.

• g (sympy expressions) – Functions on phase space.

• vars_phase (list of symbols, optional) – Phase space variables in order [x₁, p₁, x₂, p₂, …]. If None, attempt to infer from f and g.

Returns:

Poisson bracket {f, g}.

Return type:

sympy expression

symplectic.project(trajectory, plane='xy', vars_phase=None)

Project a 4D trajectory onto a 2D plane for visualization.

For 2-DOF systems, the full phase space is 4-dimensional (x₁, p₁, x₂, p₂), which cannot be directly visualized. This function extracts 2D projections suitable for plotting.

Available projection planes:

• ‘xy’ or ‘config’: Configuration space (x₁, x₂)

• ‘xp’: Phase space of first DOF (x₁, p₁)

• ‘pp’ or ‘momentum’: Momentum space (p₁, p₂)

• ‘x1p2’: Mixed coordinates (x₁, p₂)

• ‘x2p1’: Mixed coordinates (x₂, p₁)

Parameters:

• trajectory (dict) – Trajectory dictionary from hamiltonian_flow with variable name keys.

• plane (str) – Projection plane specification (default: ‘xy’). Options: ‘xy’, ‘config’, ‘xp’, ‘pp’, ‘momentum’, ‘x1p2’, ‘x2p1’

• vars_phase (list of sympy.Symbol, optional) – Variables [x1, p1, x2, p2]. If None, uses default canonical names.

Returns:

(x_data, y_data, labels) where: - x_data : ndarray, coordinates for horizontal axis - y_data : ndarray, coordinates for vertical axis - labels : tuple of str, axis label strings

Return type:

tuple

Raises:

ValueError – If the system is not 2-DOF or if plane is not recognized.

Examples

>>> traj = hamiltonian_flow(H, z0, (0, 100), vars_phase=[x1, p1, x2, p2])
>>> x, y, labels = project(traj, plane='xy')
>>> plt.plot(x, y)
>>> plt.xlabel(labels[0]); plt.ylabel(labels[1])

>>> # Phase space of first oscillator
>>> x, p, labels = project(traj, plane='xp')

Notes

Different projections reveal different aspects of the dynamics: - Configuration space shows spatial trajectories - Phase space projections show energy exchange between DOFs - Mixed projections can reveal correlations between variables

See also

visualize_poincare_section

2D section visualization

phase_portrait

1-DOF phase space visualization

symplectic.rectangle_region(center, width, height, n_points=50)

Generate a closed rectangular region in phase space.

Parameters:

• center ((float, float)) – Coordinates (x0, p0) of the rectangle centre.

• width (float) – Total width in the x direction.

• height (float) – Total height in the p direction.

• n_points (int, optional) – Number of distinct points on the perimeter (default 50). Must be at least 4.

Returns:

Ordered points forming a closed rectangle (first and last point identical). The array has shape (n_points + 1, 2).

Return type:

(N,2) ndarray

symplectic.separatrix_analysis(H, x_range, p_range, saddle_point, vars_phase=None)

Analyze separatrix structure near saddle point (1‑DOF).

Parameters:

• H (sympy expression) – Hamiltonian H(x, p).

• x_range (tuple) – Domain for visualization.

• p_range (tuple) – Domain for visualization.

• saddle_point (tuple) – Location of saddle (x_s, p_s).

• vars_phase (list, optional) – Phase space variables [x, p].

Returns:

Contains stable/unstable manifolds and energy at saddle.

Return type:

dict

symplectic.symplectic_gradient(f, vars_phase=None, numeric=False)

Compute the Hamiltonian vector field (symplectic gradient) of a function f.

For a function f on a 2n‑dimensional phase space with canonical coordinates (x₁, p₁, …, xₙ, pₙ), the symplectic gradient X_f is the unique vector field satisfying ω(X_f, ·) = df, where ω = Σ dxᵢ ∧ dpᵢ is the symplectic form. In coordinates, X_f = J⁻¹·∇f, with J⁻¹ the Poisson tensor (the inverse of the symplectic matrix).

Parameters:

• f (sympy.Expr) – Function defined on phase space.

• vars_phase (list of sympy.Symbol, optional) – Ordered list of canonical variables [x₁, p₁, x₂, p₂, …]. If None, the variables are inferred automatically from f.

• numeric (bool, default False) – If True, return a callable function that evaluates the vector field numerically at a given point. If False, return a list of sympy expressions representing each component of X_f.

Returns:

• If numeric=False (list of sympy.Expr) – Components of the Hamiltonian vector field in the same order as vars_phase, i.e., [X_f_x₁, X_f_p₁, X_f_x₂, X_f_p₂, …].

• If numeric=True (callable) – A function X(z) that takes a phase space point z (array‑like of length 2n) and returns a numpy array of the vector field components at that point.

Notes

For large‑scale evaluations (e.g., plotting on a grid), it is more efficient to create vectorized functions from the symbolic components manually:

X_components = symplectic_gradient(f, vars_phase, numeric=False)

funcs = [lambdify(vars_phase, expr, ‘numpy’) for expr in X_components]

Then evaluate on entire arrays:

results = [func(*(grid_arrays)) for func in funcs]

Examples

>>> x, p = symbols('x p', real=True)
>>> H = p**2/2 + x**2/2
>>> symplectic_gradient(H, [x, p])
[p, -x]   # dx/dt = p, dp/dt = -x

>>> # Numeric evaluation at a single point
>>> X = symplectic_gradient(H, [x, p], numeric=True)
>>> X([1.0, 0.5])
array([ 0.5, -1. ])

>>> # Vectorized evaluation for a grid (fast)
>>> X_x, X_p = symplectic_gradient(H, [x, p], numeric=False)
>>> X_x_func = lambdify((x, p), X_x, 'numpy')
>>> X_p_func = lambdify((x, p), X_p, 'numpy')
>>> q_vals = np.linspace(-2, 2, 100)
>>> p_vals = np.linspace(-2, 2, 100)
>>> Q, P = np.meshgrid(q_vals, p_vals)
>>> U = X_x_func(Q, P)   # entire array at once
>>> V = X_p_func(Q, P)

See also

poisson_bracket

{f, g} = X_f(g) = -X_g(f)

hamiltonian_flow

integrates the Hamiltonian vector field of H

SymplecticForm

the underlying symplectic structure

symplectic.visualize_phase_space_structure(H, x_range, p_range, vars_phase=None, fixed_points=None, show_separatrices=True, n_trajectories=10)

Comprehensive visualization of phase space structure for 1‑DOF.

Parameters:

• H (sympy expression) – Hamiltonian.

• x_range (tuple) – Domain ranges.

• p_range (tuple) – Domain ranges.

• vars_phase (list, optional) – Phase space variables [x, p].

• fixed_points (list, optional) – List of fixed points to analyze.

• show_separatrices (bool) – Whether to plot separatrices.

• n_trajectories (int) – Number of sample trajectories.

symplectic.visualize_poincare_section(H, z0_list, Sigma_def, vars_phase=None, tmax=100, n_returns=500, plot_vars=('x1', 'p1'), n_workers=None)

Visualize Poincaré section for multiple initial conditions (2-DOF).

Initial conditions are processed in parallel using ProcessPoolExecutor, giving a near-linear speedup with the number of available CPU cores.

Parameters:

• H (sympy.Expr) – Hamiltonian for a 2-degree-of-freedom system.

• z0_list (list of array_like) – List of initial conditions, each of length 4.

• Sigma_def (dict) – Section surface definition (see poincare_section).

• vars_phase (list of sympy.Symbol, optional) – Variables [x1, p1, x2, p2].

• tmax (float) – Maximum integration time (default: 100).

• n_returns (int) – Maximum crossings per initial condition (default: 500).

• plot_vars (tuple of str) – Two variable names to plot (default: (‘x1’, ‘p1’)).

• n_workers (int, optional) – Number of worker processes. Defaults to the number of CPU cores. Set to 1 to disable parallelism (useful for debugging).