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.

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).

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.).

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

Bases: object

Symplectic structure for an arbitrary number of degrees of freedom.

Represents the symplectic 2‑form ω = Σ dxᵢ ∧ dpᵢ (canonical) on a 2n‑dimensional phase space. The matrix representation is block‑diagonal with 2×2 blocks [[0,-1],[1,0]].

Parameters:

• n (int, optional) – Number of degrees of freedom. If not given, inferred from vars_phase.

• vars_phase (list of symbols, optional) – Ordered list of canonical variables (x1, p1, x2, p2, …).

• omega_matrix (sympy Matrix, optional) – Custom 2n×2n antisymmetric matrix. If None, canonical form is used.

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.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'))

First return map from Poincaré section points.

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)

Integrate Hamilton’s equations for any number of degrees of freedom.

Parameters:

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

• z0 (array_like) – Initial condition, length 2n.

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

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

• integrator ({'symplectic', 'verlet', 'rk45'}) – Integration method.

• n_steps (int) – Number of time steps.

Returns:

Trajectory with keys ‘t’, each variable name (e.g., ‘x1’, ‘p1’, …), and ‘energy’.

Return type:

dict

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, vars_phase=None, n_vectors=4, renorm_interval=10)

Estimate Lyapunov exponents (simplified) for 4D trajectory.

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

Monodromy matrix for a periodic orbit (2‑DOF).

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

Generate phase portrait for 1‑DOF system.

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

Poincaré section for 2‑DOF systems.

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.

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.test_coupled_oscillators()

Test 2‑DOF coupled oscillators.

symplectic.test_fixed_points_1d()

Test fixed point finding in 1‑DOF.

symplectic.test_harmonic_oscillator()

Test 1‑DOF harmonic oscillator.

symplectic.test_monodromy_simple()

Test monodromy matrix on a simple periodic orbit (2‑DOF).

symplectic.test_poincare_section()

Test Poincaré section on 2‑DOF.

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'))

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