Hamiltonian Systems

Hamiltonian mechanics provides a powerful geometric framework for describing the evolution of physical systems. Introduced by Hamilton (1834) as a reformulation of Lagrangian mechanics, Hamiltonian systems describe the dynamics of a broad range of physical phenomena—from classical mechanics and robotics to fluid dynamics and quantum mechanics.

Mathematical Foundation

A Hamiltonian system is a dynamical system evolving on a symplectic manifold \((\mathcal{M}, \omega)\) according to a smooth Hamiltonian function \(H : \mathcal{M} \to \mathbb{R}\). The symplectic form \(\omega\) is a closed, non-degenerate, differential 2-form that endows the phase space with a rich geometric structure.

For canonical Hamiltonian systems, the phase space \(\mathcal{M}\) can be modeled on the cotangent bundle \(T^*Q\) of a smooth \(n\)-dimensional configuration manifold \(Q\). In canonical coordinates \((q, p)\), where \(q \in Q\) represents the configuration and \(p \in T_q^*Q\) represents the conjugate momentum, the symplectic form takes the canonical form

which in matrix representation is given by the canonical Poisson tensor

with \(J_{2n}^\top = J_{2n}^{-1} = -J_{2n}\).

Hamilton's Equations

The evolution of a Hamiltonian system is governed by Hamilton's equations, which define the Hamiltonian vector field \(X_H = \omega^{-1} \mathrm{d}H\). For canonical systems, these reduce to the familiar first-order differential equations:

These equations describe how the configuration \(q\) and momentum \(p\) evolve in time, with the first equation relating velocity to momentum and the second expressing the force as the negative gradient of the Hamiltonian with respect to configuration.

For mechanical systems with quadratic kinetic energy structure, the Hamiltonian function decomposes into kinetic and potential energy terms:

where \(M(q)\) is the mass-inertia matrix (which may depend on the configuration \(q\) for systems with configuration-dependent inertia) and \(V(q)\) is the potential energy. The momentum is related to velocity through the Legendre transform: \(p = M(q) \dot{q}\).

Symplectic Structure and Conservation

The symplectic structure of Hamiltonian systems gives rise to fundamental conservation properties. A diffeomorphism \(f : (\mathcal{M}, \omega) \to (\mathcal{N}, \eta)\) between symplectic manifolds is a symplectomorphism if it preserves the symplectic form, i.e., \(f^*\eta = \omega\).

The Hamiltonian flow preserves the symplectic structure, ensuring that:

• Energy conservation: The Hamiltonian \(H\) is preserved along trajectories, \(\frac{\mathrm{d}H}{\mathrm{d}t} = 0\), as a direct consequence of the skew-symmetry of the symplectic form.
• Long-term stability: The symplectic structure preservation ensures that numerical integration maintains the geometric properties of the system, preventing artificial energy drift and ensuring physically consistent long-term dynamics.

Structure-Preserving Simulation in SOFAx

The SOFAx framework implements Hamiltonian systems through a modular architecture that preserves the underlying geometric structure. The design philosophy centers on Energy Fields—modular components that evaluate energies and their gradients, composing naturally to form complete Hamiltonian systems.

Computational Framework: JAX

SOFAx is built on JAX, a high-performance numerical computing library that provides three key computational advantages:

• Just-In-Time (JIT) Compilation: All SOFAx components are JIT-compilable using jax.jit or Equinox's eqx.filter_jit, enabling significant performance improvements. The entire integration loop can be compiled into optimized machine code, eliminating Python overhead and enabling aggressive compiler optimizations.

• GPU Acceleration: JAX operations automatically execute on GPU when available, and SOFAx leverages jax.lax.scan to compile entire integration loops into single GPU kernels. This enables efficient parallel execution of time-stepping algorithms, making large-scale simulations computationally tractable.

• Automatic Differentiation: Forces and constraint Jacobians are computed using JAX's automatic differentiation (jax.grad, jax.jacfwd), ensuring that derivatives are exact and preserving the geometric structure of Hamiltonian systems. This eliminates manual derivative computation and reduces the risk of implementation errors.

These computational features make SOFAx suitable for both interactive prototyping and high-performance production simulations, seamlessly scaling from CPU to GPU execution without code changes.

Energy Fields Architecture

Rather than working directly with forces and accelerations, SOFAx v3 employs Energy Fields that encapsulate the energetic structure of physical systems:

• Kinetic Energy Field: Evaluates \(T(q, p) = \frac{1}{2} p^\top M^{-1}(q) p\) where \(M(q)\) is the mass-inertia operator, which may depend on the configuration \(q\) (e.g., for systems with configuration-dependent inertia). The field provides matrix-free evaluation of kinetic energy and velocity computation through \(M^{-1}(q)\).

• Potential Energy Fields: Represent conservative forces through potential energy functions \(V(q)\). Multiple potential fields can be composed additively, enabling multi-material and multi-physics simulations.

• Dissipation Fields: Model non-conservative effects through Rayleigh dissipation functions \(D(v) = \frac{1}{2} v^\top B v\), where \(B\) is a positive semi-definite damping operator. Dissipation fields are implemented in sofax.ports.port_hamiltonian.DissipationField and used via the port-Hamiltonian extension (PortHamiltonianSystem), extending the base Hamiltonian framework to include energy-dissipating systems through the port-Hamiltonian formalism.

The Hamiltonian is composed as \(H(q, p) = T(q, p) + \sum_i V_i(q)\), giving rise to the conservative equations of motion:

For example, a complete Hamiltonian for a mechanical system with gravitational potential and elastic potential might be:

where \(m g h(q)\) represents gravitational potential energy and \(\frac{1}{2} k \|q - q_0\|^2\) represents elastic potential energy with stiffness \(k\) and rest configuration \(q_0\).

For systems with dissipation, the port-Hamiltonian extension modifies the equations of motion to:

where the dissipation term \(-\frac{\partial D}{\partial \dot{q}} = -B \dot{q}\) introduces energy dissipation into the system.

Symplectic Integration

SOFAx employs structure-preserving integrators, most notably RATTLE (an extension of the Verlet algorithm), which maintains the symplectic structure during numerical integration. RATTLE handles holonomic constraints \(\Phi(q) = 0\) through symplectic projection, ensuring that:

• Constraints are satisfied exactly at each time step
• The symplectic structure is preserved
• Energy conservation is maintained for conservative systems
• Long-term stability is guaranteed

The integration scheme alternates between momentum updates (kicks) and position updates (drifts), with constraint projections applied at appropriate stages to maintain both the symplectic structure and constraint satisfaction.

The RATTLE Algorithm

RATTLE (Rychaert–Andersen–Teller–Loupis–Ermak) is a second-order symplectic integrator that extends the Verlet algorithm to constrained Hamiltonian systems. RATTLE is itself an extension of the SHAKE algorithm, adding velocity constraint enforcement to maintain symplecticity. Given a time step \(\Delta t\) and holonomic constraints \(\Phi(q) = 0\) with Jacobian \(G(q) = \partial \Phi / \partial q\), RATTLE advances the system from \((q^n, p^n)\) to \((q^{n+1}, p^{n+1})\) through five stages:

1. Pre-kick: Update momentum by half a step using the conservative forces:

Note: For conservative systems, the force \(-\frac{\partial H}{\partial q}\) depends only on the configuration \(q^n\), not on the momentum \(p^n\). The momentum parameter is included in the API for consistency and future support of dissipative systems where forces may depend on both \(q\) and \(p\).

1. Drift: Update position using the intermediate momentum:

1. Position projection: Project \(q^*\) onto the constraint manifold \(\Phi(q) = 0\) to obtain \(q^{n+1}\), ensuring exact constraint satisfaction. This projection is performed using Newton-Raphson iteration with pseudo-inverse on the constraint violation, with the mass-weighted metric preserving the symplectic structure.

2. Post-kick: Complete the momentum update using the new position:

For conservative systems, the force \(-\frac{\partial H}{\partial q}\) depends only on the configuration \(q\), not on the momentum \(p\). For dissipative systems (handled by the port-Hamiltonian extension), the force may also depend on \(p\).

1. Velocity projection: Project \(p^*\) to satisfy the velocity constraint \(G(q^{n+1}) \cdot M^{-1} p = 0\), yielding the final momentum \(p^{n+1}\).

The algorithm's symplecticity stems from two key features: (1) the symmetric structure of the kicks (half-step before and after the drift), which mirrors the structure of the continuous Hamiltonian flow, and (2) symplectic projections that preserve the canonical symplectic form \(\omega\) during constraint enforcement. This construction ensures that the discrete flow approximates the continuous Hamiltonian flow while maintaining the geometric structure of the phase space, guaranteeing long-term stability and energy conservation for conservative systems.

Constraint Projection Methods

The projection steps in RATTLE (steps 3 and 5) are critical for maintaining both constraint satisfaction and symplecticity. In SOFAx, these projections are handled through specialized methods in constraint field classes, which define the constraint function \(\Phi(q)\) and its Jacobian \(G(q) = \partial \Phi / \partial q\).

Each constraint field implements two projection methods that work together to enforce constraints at both the position and velocity levels:

Position Projection (project_position): This method projects an unconstrained configuration \(q^*\) onto the constraint manifold \(\Phi(q) = 0\) to obtain \(q^{n+1}\). The implementation strategy depends on the constraint type:

• Fixed-point constraints: Direct assignment—the constrained degrees of freedom are set to their prescribed values. This is a non-iterative, exact operation.

• General holonomic constraints: Iterative Newton-Raphson correction with pseudo-inverse. The method minimizes the constraint violation \(\|\Phi(q)\|\) while preserving the symplectic structure through a mass-weighted metric, ensuring that the projection respects the geometric structure of the phase space.

Velocity Projection (project_velocity): This method projects momentum \(p^*\) onto the tangent space of the constraint manifold, ensuring that the velocity \(v = M^{-1} p\) satisfies the velocity constraint \(G(q) \cdot v = 0\). Like position projection, the strategy depends on the constraint type:

• Fixed-point constraints: Direct (non-iterative) projection—the momentum at constrained degrees of freedom is set to zero, corresponding to zero velocity at those points.

• General constraints: Iterative symplectic projection that adjusts the momentum to satisfy \(G(q) \cdot M^{-1} p = 0\) while preserving the canonical symplectic structure.

The modular design of constraint fields allows each constraint type to implement its own efficient projection strategy, from simple fixed-point assignments to complex iterative projections for volumetric or geometric constraints. This flexibility enables SOFAx to handle diverse constraint types—from simple boundary conditions to complex volume-preserving constraints—while maintaining the fundamental geometric properties of Hamiltonian systems.

By embedding the symplectic structure directly into both the integration scheme and the constraint projection methods, SOFAx guarantees that numerical solutions preserve the fundamental conservation laws and geometric properties of the underlying physical system. This structure-preserving approach enables accurate and stable long-term simulations of Hamiltonian dynamics, making it particularly well-suited for applications requiring energy conservation and long-time stability.