[Paper Review] On converting any one-step method to a variational integrator of the same order

George W. Patrick, Raymond J. Spiteri+2 more|ArXiv.org|Mar 5, 2009

Numerical methods for differential equations15 references3 citations

TL;DR

This paper presents a general algorithm to convert any one-step numerical integrator into a variational integrator of the same order for constrained mechanical systems. By leveraging the discrete variational principle and constructing a symplectic layer using the standard integrator's output, the method ensures structure-preserving properties like symplecticity and momentum conservation while maintaining the original method's order of accuracy.

ABSTRACT

In the formalism of constrained mechanics, such as that which underlies the SHAKE and RATTLE methods of molecular dynamics, we present an algorithm to convert any one-step integration method to a variational integrator of the same order. The one-step method is arbitrary, and the conversion can be automated, resulting in a powerful and flexible approach to the generation of novel variational integrators with arbitrary order.

Motivation & Objective

• To develop a systematic method for transforming arbitrary one-step numerical integrators into variational integrators of the same order for constrained mechanical systems.
• To preserve key geometric properties such as symplecticity and momentum conservation in the resulting integrator.
• To enable automated generation of novel high-order variational integrators from existing standard integrators.
• To ensure the conversion process maintains the same order of accuracy as the original one-step method.
• To provide a framework applicable to systems with holonomic constraints, including molecular dynamics and rigid body dynamics.

Proposed method

• The method uses a standard one-step integrator as a 'standard layer' to generate curve segments and approximate action values for discrete Lagrangian construction.
• A symplectic layer is constructed by defining bias parameters α⁻ and α⁺ to control the discrete time evolution and ensure consistency with the variational principle.
• The discrete Euler–Lagrange equations are derived using approximations of the discrete Lagrangian and its derivatives, with key approximations based on the standard layer's output.
• The algorithm solves a system of linear equations involving constraint forces and multipliers, using matrix forms such as [I + D²(θᵀg), Dgᵀ; Dg, 0] to enforce consistency.
• Iterative solution of the nonlinear system is performed using approximations like (4.8), which are computationally equivalent to solving the full linear system.
• Constraint forces are tracked via multipliers ν₁, ν₂, which are non-unique but whose physical interpretations (e.g., force of constraint) are preserved in products like ν₂ᵀDg(q).

Experimental results

Research questions

• RQ1 Can any one-step numerical integrator be systematically converted into a variational integrator of the same order?
• RQ2 How can symplecticity and momentum conservation be preserved during the conversion of a standard integrator?
• RQ3 What is the role of bias parameters α⁻ and α⁺ in constructing a consistent discrete variational integrator?
• RQ4 How can the discrete variational principle be enforced when the original integrator is not inherently symplectic?
• RQ5 What conditions ensure that the resulting integrator maintains the same order of accuracy as the original one-step method?

Key findings

• The proposed algorithm successfully converts any one-step integrator into a variational integrator of the same order, preserving symplecticity and momentum conservation.
• The method is fully automated and applicable to arbitrary one-step methods, enabling the generation of novel high-order variational integrators.
• The discrete Lagrangian is constructed using the standard layer’s output, ensuring consistency with the classical action up to the desired order.
• The solution of the discrete variational equations relies on solving a linear system approximated via (4.8), which is computationally equivalent to the exact system.
• Constraint forces are encoded in non-unique multipliers, but their physical meaning is preserved in products such as ν₂ᵀDg(q), which represent the actual force of constraint.
• The iterative solution process converges to the correct state variables q₂ and v₂, with convergence monitored via physical quantities rather than the multipliers themselves.