Structure-preserving Integrators Research Articles

Hamiltonian formulation of X-point collapse in an extended magnetohydrodynamics framework

The study of X-point collapse in magnetic reconnection has witnessed extensive research in the context of space and laboratory plasmas. In this paper, a recently derived mathematical formulation of X-point collapse applicable in the regime of extended magnetohydrodynamics is shown to possess a noncanonical Hamiltonian structure composed of five dynamical variables inherited from its parent model. The Hamiltonian and the noncanonical Poisson brackets are both derived, and the latter is shown to obey the requisite properties of antisymmetry and the Jacobi identity (an explicit proof of the latter is provided). In addition, the governing equations for the Casimir invariants are presented, and one such solution is furnished. The above features can be harnessed and expanded in future work, such as developing structure-preserving integrators for this dynamical system.

Shadow Hamiltonians of structure-preserving integrators for Nambu mechanics

Abstract Symplectic integrators are widely implemented numerical integrators for Hamiltonian mechanics, which preserve the Hamiltonian structure (symplecticity) of the system. Although the symplectic integrator does not conserve the energy of the system, it is well known that there exists a conserving modified Hamiltonian, called the shadow Hamiltonian. For the Nambu mechanics, which is a kind of generalized Hamiltonian mechanics, we can also construct structure-preserving integrators by the same procedure used to construct the symplectic integrators. In the structure-preserving integrator, however, the existence of shadow Hamiltonians is nontrivial. This is because the Nambu mechanics is driven by multiple Hamiltonians and it is nontrivial whether the time evolution by the integrator can be cast into the Nambu mechanical time evolution driven by multiple shadow Hamiltonians. In this paper we present a general procedure to calculate the shadow Hamiltonians of structure-preserving integrators for Nambu mechanics, and give an example where the shadow Hamiltonians exist. This is the first attempt to determine the concrete forms of the shadow Hamiltonians for a Nambu mechanical system. We show that the fundamental identity, which corresponds to the Jacobi identity in Hamiltonian mechanics, plays an important role in calculating the shadow Hamiltonians using the Baker–Campbell–Hausdorff formula. It turns out that the resulting shadow Hamiltonians have indefinite forms depending on how the fundamental identities are used. This is not a technical artifact, because the exact shadow Hamiltonians obtained independently have the same indefiniteness.

Principled and interpretable alignability testing and integration of single-cell data

Single-cell data integration can provide a comprehensive molecular view of cells, and many algorithms have been developed to remove unwanted technical or biological variations and integrate heterogeneous single-cell datasets. Despite their wide usage, existing methods suffer from several fundamental limitations. In particular, we lack a rigorous statistical test for whether two high-dimensional single-cell datasets are alignable (and therefore should even be aligned). Moreover, popular methods can substantially distort the data during alignment, making the aligned data and downstream analysis difficult to interpret. To overcome these limitations, we present a spectral manifold alignment and inference (SMAI) framework, which enables principled and interpretable alignability testing and structure-preserving integration of single-cell data with the same type of features. SMAI provides a statistical test to robustly assess the alignability between datasets to avoid misleading inference and is justified by high-dimensional statistical theory. On a diverse range of real and simulated benchmark datasets, it outperforms commonly used alignment methods. Moreover, we show that SMAI improves various downstream analyses such as identification of differentially expressed genes and imputation of single-cell spatial transcriptomics, providing further biological insights. SMAI's interpretability also enables quantification and a deeper understanding of the sources of technical confounders in single-cell data.

Structure-preserving integrators based on a new variational principle for constrained mechanical systems

A new variational principle for mechanical systems subject to holonomic constraints is presented. The newly proposed GGL principle is closely related to the often used Gear-Gupta-Leimkuhler (GGL) stabilization of the differential–algebraic equations governing the motion of constrained mechanical systems. The GGL variational principle relies on an extension of the Livens principle (sometimes also referred to as Hamilton–Pontryagin principle) to mechanical systems subject to holonomic constraints. In contrast to the original GGL stabilization, the new approach facilitates the design of structure-preserving integrators. In particular, new variational integrators are presented, which result from the direct discretization of the GGL variational principle. These variational integrators are symplectic and conserve momentum maps in the case of systems with symmetry. In addition to that, a new energy–momentum scheme is developed, which results from the discretization of the Euler–Lagrange equations pertaining to the GGL variational principle. The numerical properties of the newly devised schemes are investigated in representative examples of constrained mechanical systems.

A structure-preserving integrator for incompressible finite elastodynamics based on a grad-div stabilized mixed formulation with particular emphasis on stretch-based material models

We present a structure-preserving scheme based on a recently-proposed mixed formulation for incompressible hyperelasticity formulated in principal stretches. Although there exist several different Hamiltonians introduced for quasi-incompressible elastodynamics based on different multifield variational formulations, there is not much study on the fully incompressible materials in the literature. The adopted mixed formulation can be viewed as a finite-strain generalization of Herrmann variational formulation, and it naturally provides a new Hamiltonian for fully incompressible elastodynamics. Invoking the discrete gradient and scaled mid-point formulas, we are able to design fully-discrete schemes that preserve the Hamiltonian and momenta. Our analysis and numerical evidence also reveal that the scaled mid-point formula is non-robust numerically. The generalized Taylor–Hood element based on the spline technology conveniently provides a higher-order, robust, and inf-sup stable spatial discretization option for finite strain analysis. To enhance the element performance in volume conservation, the grad-div stabilization, a technique initially developed in computational fluid dynamics, is introduced here for elastodynamics. It is shown that the stabilization term does not impose additional restrictions for the algorithmic stress to respect the invariants, leading to an energy-decaying and momentum-conserving fully discrete scheme. A set of numerical examples is provided to justify the claimed properties. The grad-div stabilization is found to enhance the discrete mass conservation effectively. Furthermore, in contrast to conventional algorithms based on Cardano’s formula and perturbation techniques, the spectral decomposition algorithm developed by Scherzinger and Dohrmann is robust and accurate to ensure the discrete conservation laws and is thus recommended for stretch-based material modeling.

On the design of non-singular, energy-momentum consistent integrators for nonlinear dynamics using energy splitting and perturbation techniques

This work proposes a suite of numerical techniques to facilitate the design of structure-preserving integrators for nonlinear dynamics, with particular emphasis on many-body dynamics. The celebrated LaBudde-Greenspan integrator and various energy-momentum schemes adopt a difference quotient formula in their algorithmic force definitions, which suffers from numerical instability as the denominator gets close to zero. There is a need to develop structure-preserving integrators without invoking the quotient formula. In this work, the potential energy of a Hamiltonian system is split into two parts, and specially developed quadrature rules are applied separately to them. The resulting integrators can be regarded as classical ones perturbed with first- or second-order terms, and the energy split guarantees the dissipative nature in the numerical residual. In the meantime, the conservation of invariants is respected in the design. A complete analysis of the proposed integrators is given, with representative numerical examples provided to demonstrate their performance. They can be used either independently as energy-decaying and momentum-conserving schemes for nonlinear problems or as an alternate option with a conserving integrator, such as the LaBudde-Greenspan integrator, when the numerical instability in the difference quotient is detected.

The GGL variational principle for constrained mechanical systems

We present an extension of the Livens variational principle (sometimes also referred to as Hamilton-Pontryagin principle) to mechanical systems subject to holonomic constraints. The newly proposed principle embodies an index reduction in the spirit of the often-applied GGL stabilization and thus may be termed “GGL principle”. The Euler-Lagrange equations of the GGL principle assume the form of differential-algebraic equations (DAEs) with differentiation index two. In contrast to the original GGL-DAEs, the present formulation fits into the Hamiltonian framework of mechanics. Therefore, the GGL principle facilitates the design of structure-preserving integrators. In particular, it offers the possibility to construct variational integrators. This is illustrated with the development of a new first-order scheme which is symplectic by design. The numerical properties of the newly devised scheme are investigated for representative examples of constrained mechanical systems.

NySALT: Nyström-type inference-based schemes adaptive to large time-stepping

Large time-stepping is important for efficient long-time simulations of deterministic and stochastic Hamiltonian dynamical systems. Conventional structure-preserving integrators, while being successful for generic systems, have limited tolerance to time step size due to stability and accuracy constraints. We propose to use data to innovate classical integrators so that they can be adaptive to large time-stepping and are tailored to each specific system. In particular, we introduce NySALT, Nyström-type inference-based schemes adaptive to large time-stepping. The NySALT has optimal parameters for each time step learnt from data by minimizing the one-step prediction error. Thus, it is tailored for each time step size and the specific system to achieve optimal performance and tolerate large time-stepping in an adaptive fashion. We prove and numerically verify the convergence of the estimators as data size increases. Furthermore, analysis and numerical tests on the deterministic and stochastic Fermi-Pasta-Ulam (FPU) models show that NySALT enlarges the maximal admissible step size of linear stability, and quadruples the time step size of the Störmer–Verlet and the BAOAB when maintaining similar levels of accuracy.

Commutator-free Lie group methods with minimum storage requirements and reuse of exponentials

A new format for commutator-free Lie group methods is proposed based on explicit classical Runge-Kutta schemes. In this format exponentials are reused at every stage and the storage is required only for two quantities: the right hand side of the differential equation evaluated at a given Runge-Kutta stage and the function value updated at the same stage. The next stage of the scheme is able to overwrite these values. The result is proven for a 3-stage third order method and a conjecture for higher order methods is formulated. Five numerical examples are provided in support of the conjecture. This new class of structure-preserving integrators has a wide variety of applications for numerically solving differential equations on manifolds.

On dissipative symplectic integration with applications to gradient-based optimization

Recently, continuous-time dynamical systems have proved useful in providing conceptual and quantitative insights into gradient-based optimization, widely used in modern machine learning and statistics. An important question that arises in this line of work is how to discretize the system in such a way that its stability and rates of convergence are preserved. In this paper we propose a geometric framework in which such discretizations can be realized systematically, enabling the derivation of ‘rate-matching’ algorithms without the need for a discrete convergence analysis. More specifically, we show that a generalization of symplectic integrators to non-conservative and in particular dissipative Hamiltonian systems is able to preserve rates of convergence up to a controlled error. Moreover, such methods preserve a shadow Hamiltonian despite the absence of a conservation law, extending key results of symplectic integrators to non-conservative cases. Our arguments rely on a combination of backward error analysis with fundamental results from symplectic geometry. We stress that although the original motivation for this work was the application to optimization, where dissipative systems play a natural role, they are fully general and not only provide a differential geometric framework for dissipative Hamiltonian systems but also substantially extend the theory of structure-preserving integration.

Asymptotic analysis of a structure-preserving integrator for damped Hamiltonian systems

The present work deals with the numerical long-time integration of damped Hamiltonian systems. The method that we analyze combines a specific Strang splitting, that separates linear dissipative effects from conservative ones, with an energy-preserving averaged vector field (AVF) integrator for the Hamiltonian subproblem. This construction faithfully reproduces the energy-dissipation structure of the continuous model, its equilibrium points and its natural Lyapunov function. As a consequence of these structural similarities, both the convergence to equilibrium and, more interestingly, the energy decay rate of the continuous dynamical system are recovered at a discrete level. The possibility of replacing the implicit AVF integrator by an explicit Stormer-Verlet one is also discussed, while numerical experiments illustrate and support the theoretical findings.

Structure-preserving integrators for dissipative systems based on reversible- irreversible splitting.

We study the optimal design of numerical integrators for dissipative systems, for which there exists an underlying thermodynamic structure known as GENERIC (general equation for the nonequilibrium reversible–irreversible coupling). We present a frame-work to construct structure-preserving integrators by splitting the system into reversible and irreversible dynamics. The reversible part, which is often degenerate and reduces to a Hamiltonian form on its symplectic leaves, is solved by using a symplectic method (e.g. Verlet) with degenerate variables being left unchanged, for which an associated modified Hamiltonian (and subsequently a modified energy) in the form of a series expansion can be obtained by using backward error analysis. The modified energy is then used to construct a modified friction matrix associated with the irreversible part in such a way that a modified degeneracy condition is satisfied. The modified irreversible dynamics can be further solved by an explicit midpoint method if not exactly solvable. Our findings are verified by various numerical experiments, demonstrating the superiority of structure-preserving integrators over alternative schemes in terms of not only the accuracy control of both energy conservation and entropy production but also the preservation of the conformal symplectic structure in the case of linearly damped systems.

A review of some geometric integrators

Some of the most important geometric integrators for both ordinary and partial differential equations are reviewed and illustrated with examples in mechanics. The class of Hamiltonian differential systems is recalled and its symplectic structure is highlighted. The associated natural geometric integrators, known as symplectic integrators, are then presented. In particular, their ability to numerically reproduce first integrals with a bounded error over a long time interval is shown. The extension to partial differential Hamiltonian systems and to multisymplectic integrators is presented afterwards. Next, the class of Lagrangian systems is described. It is highlighted that the variational structure carries both the dynamics (Euler–Lagrange equations) and the conservation laws (Nœther’s theorem). Integrators preserving the variational structure are constructed by mimicking the calculus of variation at the discrete level. We show that this approach leads to numerical schemes which preserve exactly the energy of the system. After that, the Lie group of local symmetries of partial differential equations is recalled. A construction of Lie-symmetry-preserving numerical scheme is then exposed. This is done via the moving frame method. Applications to Burgers equation are shown. The last part is devoted to the Discrete Exterior Calculus, which is a structure-preserving integrator based on differential geometry and exterior calculus. The efficiency of the approach is demonstrated on fluid flow problems with a passive scalar advection.

Discrete-time port-Hamiltonian systems based on Gauss-Legendre collocation

We introduce a family of discrete-time lossless input-state-output port-Hamiltonian systems based on numerical time integration with symplectic collocation schemes. For systems with non-zero input, symplecticity extends to the conservation of a discrete energy balance, based on which a discrete-time Dirac structure is defined. Using Gauss-Legendre collocation, the corresponding quadrature formula allows to quantify the discretization error for the supplied energy. On a linear example, backward error analysis and numerical experiments are performed in order to illustrate the accuracy of the resulting structure-preserving integration schemes.

Structure-preserving time integration of non-isothermal finite viscoelastic continua related to variational formulations of continuum dynamics

A structure-preserving time integration scheme is a modified standard time stepping scheme, which fulfills the balance laws of a dynamical problem also in a discrete setting. Hence, such a scheme represents a space-time discretization of a continuum with all of its central properties. This paper deals with a special second-order accurate structure-preserving time integration of thermo-mechanical couplings in a moving continuum, satisfying all discrete thermodynamical balance laws in contrast to the underlying second-order accurate midpoint rule. More accurately, we consider isotropic heat conduction and the Gough–Joule effect in moving continua. The considered time evolution is described by five differential equations with respect to five independent variables, including the entropy field of the continuum. The bodies in the numerical examples are subjected to Dirichlet- or Neumann boundary conditions in the displacement as well as the temperature field. Although this structure-preserving time integrator is based on a time discrete spatially weak finite element formulation, it fulfills the same balance laws as the underlying five differential equations. Namely, in addition to the balances of entropy, linear momentum and angular momentum, also the balances of total energy and Lyapunov’s function are fulfilled in the discrete case.

On the integration of singularity-free representations of $$\varvec{SO(3)}$$ S O ( 3 ) for direct optimal control

In this paper we analyze the performance of different combinations of: (1) parameterization of the rotational degrees of freedom (DOF) of multibody systems, and (2) choice of the integration scheme, in the context of direct optimal control discretized according to the direct multiple-shooting method. The considered representations include quaternions and Direction Cosine Matrices, both having the peculiarity of being non-singular and requiring more than three parameters to describe an element of the Special Orthogonal group $$\textit{SO}(3)$$ . These representations yield invariants in the dynamics of the system, i.e., algebraic conditions which have to be satisfied in order for the model to be representative of physical reality. The investigated integration schemes include the classical explicit Runge–Kutta method, its stabilized version based on Baumgarte’s technique, which tends to reduce the drift from the underlying manifold, and a structure-preserving alternative, namely the Runge–Kutta Munthe-Kaas method, which preserves the invariants by construction. The performances of the combined choice of representation and integrator are assessed by solving thousands of planning tasks for a nonholonomic, underactuated cart-pendulum system, where the pendulum can experience arbitrarily large 3D rotations. The aspects analyzed include success rate, average number of iterations and CPU time to convergence, and quality of the solution. The results reveal how structure-preserving integrators are the only choice for lower accuracies, whereas higher-order, non-stabilized standard integrators seem to be the computationally most competitive solution when higher levels of accuracy are pursued. Overall, the quaternion-based representation is the most efficient in terms of both iterations and CPU time to convergence, albeit at the cost of lower success rates and increased probability of being trapped by higher local minima.

Energy-Preserving Algorithms for the Benjamin Equation

This paper presents several hybrid algorithms to preserve the global energy of the Benjamin equation. The Benjamin equation is a non-local partial differential equation involving the Hilbert transform. For this sake, quite few structure-preserving integrators have been proposed so far. Our schemes are derived based on an extended multi-symplectic Hamiltonian system of the Benjamin equation by using Fourier pseudospectral method, finite element method and wavelet collocation method spatially coupled with the AVF method temporally. The local and global properties of the proposed schemes are studied. Numerical experiments are presented to demonstrate the conservative properties of the proposed numerical methods and study the evolutions of the numerical solutions of solitary waves and wave breaking.

Symmetric parareal algorithms for Hamiltonian systems

The parareal in time algorithm allows for efficient parallel numerical simulations of time-dependent problems. It is based on a decomposition of the time interval into subintervals, and on a predictor-corrector strategy, where the propagations over each subinterval for the corrector stage are concurrently performed on the different processors that are available. In this article, we are concerned with the long time integration of Hamiltonian systems. Geometric, structure-preserving integrators are preferably employed for such systems because they show interesting numerical properties, in particular excellent preservation of the total energy of the system. Using a symmetrization procedure and/or a (possibly also symmetric) projection step, we introduce here several variants of the original plain parareal in time algorithm [L. Baffico, et al. Phys. Rev. E 66 (2002) 057701; G. Bal and Y. Maday, A parareal time discretization for nonlinear PDE’s with application to the pricing of an American put, in Recent developments in domain decomposition methods, Lect. Notes Comput. Sci. Eng. 23 (2002) 189–202; J.-L. Lions, Y. Maday and G. Turinici, C. R. Acad. Sci. Paris, Série I 332 (2001) 661–668.] that are better adapted to the Hamiltonian context. These variants are compatible with the geometric structure of the exact dynamics, and are easy to implement. Numerical tests on several model systems illustrate the remarkable properties of the proposed parareal integrators over long integration times. Some formal elements of understanding are also provided.

Natural Coordinates in the Optimal Control of Multibody Systems

The formulation of multibody dynamics in terms of natural coordinates (NCs) leads to equations of motion in the form of differential-algebraic equations (DAEs). A characteristic feature of the natural coordinates approach is a constant mass matrix. The DAEs make possible (i) the systematic assembly of open-loop and closed-loop multibody systems, (ii) the design of state-of-the-art structure-preserving integrators such as energy-momentum or symplectic-momentum schemes, and (iii) the direct link to nonlinear finite element methods. However, the use of NCs in the optimal control of multibody systems presents two major challenges. First, the consistent application of actuating joint-forces becomes an issue since conjugate joint-coordinates are not directly available. Second, numerical methods for optimal control with index-3 DAEs are still in their infancy. The talk will address the two aforementioned issues. In particular, a new energy-momentum consistent method for the optimal control of multibody systems in terms of NCs will be presented.

Symplectic structure-preserving integrators for the two-dimensional Gross–Pitaevskii equation for BEC

Symplectic integrators have been developed for solving the two-dimensional Gross–Pitaevskii equation. The equation is transformed into a Hamiltonian form with symplectic structure. Then, symplectic integrators, including the midpoint rule, and a splitting symplectic scheme are developed for treating this equation. It is shown that the proposed codes fulfill the discrete charge conservation law. Furthermore, the global error of the numerical solution is theoretically estimated. The theoretical analysis is supported by some numerical simulations.