Discrete variational calculus for double-bracket dissipation

1. Problem Statement

The paper addresses the challenge of numerically simulating mechanical systems subject to double-bracket dissipation. These are dissipative systems where the energy decreases along the trajectory, but specific geometric structures—specifically the coadjoint orbits of the underlying Lie algebra—must remain invariant.

• Context: Standard geometric integrators (like symplectic integrators) excel at preserving geometric properties (symplecticity, momentum maps) for conservative systems. However, adding dissipation typically destroys these invariants.
• Specific Challenge: In systems like satellites with dampers, geophysical fluids, and plasma physics, the dynamics evolve on the dual of a Lie algebra ($\mathfrak{g}^*$). While energy is dissipated, the system must strictly remain on a specific coadjoint orbit (a level set of Casimir invariants).
• Limitation of Existing Methods: General-purpose integrators (e.g., Runge-Kutta) often fail to preserve the coadjoint orbit exactly, leading to drift in Casimir invariants and qualitatively incorrect long-term behavior, even if they are high-order accurate.

2. Methodology

The authors develop a Discrete Variational Integrator specifically tailored for forced Euler–Poincaré and Lie–Poisson systems with double-bracket dissipation.

A. Theoretical Foundation

1. Double-Bracket Dissipation: The paper utilizes the framework established by Bloch et al. (1996), where the dissipative force is defined via a map $\phi: \mathfrak{g} \to \mathfrak{g}$. The continuous equations of motion are:
   $$\dot{\mu} = \text{ad}^*_{\frac{\delta H}{\delta \mu}} \mu + \text{ad}^*_{\phi(\frac{\delta H}{\delta \mu})} \mu$$
   This formulation ensures that while the Hamiltonian $H$ decreases ($\dot{H} \leq 0$), the Casimir functions $C$ are conserved ($\dot{C} = 0$), keeping the trajectory on the coadjoint orbit.

2. Discrete Variational Principle:

   • The authors define a discrete Lagrangian $L_d: G \to \mathbb{R}$ on the Lie group $G$.
   • They introduce a discrete generalized variational principle with specific variations $\delta w_k$ that incorporate the dissipative map $\phi_d$.
   • The discrete action is minimized subject to variations of the form:
     $$\delta w_k = -\eta_k w_k + w_k (\text{Ad}_{\phi_d(w_k)} \eta_{k+1})$$
     where $\eta_k$ are variations in the Lie algebra vanishing at boundaries.

3. Retraction Maps:
   To construct practical algorithms, the authors use retraction maps $\tau: \mathfrak{g} \to G$ (analytic local diffeomorphisms mapping the Lie algebra to the Lie group, e.g., the Exponential map or the Cayley map).

   • The discrete momentum $M_k$ is updated via:
     $$M_{k+1} = \text{Ad}^*_{w_k \phi_d(w_k)} M_k$$
   • This update rule is derived by discretizing the continuous variational principle using the retraction map to relate group increments ($w_k$) to Lie algebra velocities ($\xi_k$).

B. Algorithm Construction

The proposed Discrete Double-Bracket (DDB) method proceeds as follows:

1. Initialize discrete momentum $M_0 = \mu(t_0)$.
2. For each time step $k$:
   • Solve for the Lie algebra element $\xi_k$ such that $M_k = (d\tau^{-1}_{h\xi_k})^* \frac{\partial l}{\partial \xi}(\xi_k)$.
   • Compute the group increment $w_k = \tau(h\xi_k)$.
   • Update momentum: $M_{k+1} = \text{Ad}^*_{w_k \phi_d(w_k)} M_k$.
3. The map $\phi_d$ is constructed from the continuous $\phi$ using the retraction: $\phi_d(w_k) = \tau(h\phi(\xi_k))$.

3. Key Contributions

1. Exact Preservation of Coadjoint Orbits: The primary contribution is a geometric integrator that exactly preserves coadjoint orbits in the presence of double-bracket dissipation. Unlike standard methods where Casimir invariants drift, this method ensures the numerical solution stays strictly on the manifold defined by the initial conditions.
2. Variational Derivation: The method is derived from a discrete variational principle, ensuring it inherits the stability and structure-preserving properties of variational integrators.
3. Flexibility via Retraction Maps: The framework allows for the use of different retraction maps (Exponential, Cayley, etc.), enabling the construction of integrators with different computational costs and orders of accuracy (e.g., second-order symmetric schemes).
4. Application to Rigid Body: The theory is explicitly applied to the dissipative rigid body, providing a concrete matrix formulation for $SO(3)$.

4. Results and Numerical Simulations

The authors validate the method using the dissipative rigid body as a test case, comparing the DDB method (using both Exponential and Cayley retractions) against:

• Runge-Kutta 4th Order (RK4)
• 3-stage Lobatto IIIC
• DOP853 (8th order, used as a reference)

Key Findings:

• Casimir Preservation: In all simulations, the DDB methods preserved the Casimir invariant (the magnitude of angular momentum $\|M\|$) to machine precision. In contrast, RK4 and Lobatto IIIC showed significant drift in the Casimir, causing the trajectory to leave the coadjoint sphere.
• Asymptotic Behavior:
  • Scenario A (Isolated Equilibria): For a rigid body with distinct moments of inertia, the system converges to rotation about the axis of maximum inertia (minimum energy). The DDB methods converged to the correct limit point, while non-geometric methods drifted due to Casimir errors.
  • Scenario B (Great Circle Equilibria): For a symmetric body (two equal moments of inertia), the system converges to a great circle of equilibria. The DDB methods correctly converged to a point on this circle, whereas other methods drifted off the manifold.
• Efficiency: The DDB methods (specifically with Cayley retraction) were computationally efficient, often faster than the high-order Lobatto IIIC method while providing superior qualitative accuracy.
• Convergence: The methods demonstrated convergence to the reference solution, with the error in the final momentum being significantly lower for DDB methods in long-time simulations due to the preservation of the underlying geometry.

5. Significance

This work bridges a critical gap in geometric numerical integration:

• Physical Fidelity: For systems like satellites with internal dampers or plasma confinement, the conservation of angular momentum (or the equivalent Casimir) is physically mandatory. Standard integrators violate this, leading to unphysical results over long simulation times.
• Robustness: The proposed integrators are robust for long-time integration, correctly capturing the asymptotic stability of the system without artificial drift.
• Generalizability: The framework is not limited to rigid bodies. The authors note that this approach can be extended to Lie algebroids and ideal fluids, offering a pathway to simulate complex fluid dynamics with dissipation while preserving the infinite-dimensional Casimir invariants of Euler's equations.

In conclusion, the paper presents a rigorous, variational-based geometric integrator that successfully handles double-bracket dissipation, ensuring that numerical simulations of dissipative mechanical systems remain geometrically faithful to the continuous physics.