Symplectic Integrators Research Articles

Computational efficiency of symplectic integration schemes: application to multidimensional disordered Klein–Gordon lattices

We implement several symplectic integrators, which are based on two part splitting, for studying the chaotic behavior of one- and two-dimensional disordered Klein-Gordon lattices with many degrees of freedom and investigate their numerical performance. For this purpose, we perform extensive numerical simulations by considering many different initial energy excitations and following the evolution of the created wave packets in the various dynamical regimes exhibited by these models. We compare the efficiency of the considered integrators by checking their ability to correctly reproduce several features of the wave packets propagation, like the characteristics of the created energy distribution and the time evolution of the systems' maximum Lyapunov exponent estimator. Among the tested integrators the fourth order $ABA864$ scheme \cite{BCFLMM13} showed the best performance as it needed the least CPU time for capturing the correct dynamical behavior of all considered cases when a moderate accuracy in conserving the systems' total energy value was required. Among the higher order schemes used to achieve a better accuracy in the energy conservation, the sixth order scheme $s11ABA82\_6$ exhibited the best performance.

A symplectic procedure for two-dimensional coupled elastic wave equations using radial basis functions interpolation

This study proposes a meshless symplectic algorithm for two-dimensional coupled elastic waves, which satisfy the equations of linear isotropic elasticity on a rectangular domain with appropriate initial and boundary conditions. This method is based on a combination of radial basis functions (RBFs) collocation and geometric integrators. The method involves two steps, first the spatial discretization by using RBFs interpolation transforms the seismic waves into a finite-dimensional Hamiltonian system, then evolves the semi-discretized system by using symplectic integrators. Furthermore, the conservations of energy and symplectic structure are investigated. The truncation error and global error of the method are estimated. Numerical experiments show that the scheme is efficient, easy to implement with the scattered knots and possesses a good behavior in the long time integration.

Symplectic Radial Basis Approximation of Multi-variate Hamiltonian PDEs

A meshless symplectic algorithm for multi-dimensional linear Hamiltonian partial differential equations with radial basis interpolation has been suggested when the solving domain is $${\mathbb {R}}^d$$ . In this study, for sake of application, the energy error estimate of the proposed symplectic scheme is made when the targeted domain is bounded. The paper also extends the algorithm to solve the nonlinear two-dimensional Klein–Gordon equation. As numerical experiments shown, the algorithm preserves the discretization energy, and it is efficient to deal with the nonlinear case.

Explicit symplectic integrator for particle tracking in s -dependent static electric and magnetic fields with curved reference trajectory

We describe a method for symplectic tracking of charged particles through static electric and magnetic fields. The method can be applied to cases where the fields have a dependence on longitudinal as well as transverse position, and where the reference trajectory may have non-zero curvature. Application of the method requires analytical expressions for the scalar and vector potentials: we show how suitable expressions, in the form of series analogous to multipole expansions, can be constructed from numerical field data, allowing the method to be used in cases where only numerical field data are available.

Systematic corrections to the Thomas–Fermi approximation without a gradient expansion

We improve on the Thomas–Fermi approximation for the single-particle density of fermions by introducing inhomogeneity corrections. Rather than invoking a gradient expansion, we relate the density to the unitary evolution operator for the given effective potential energy and approximate this operator by a Suzuki–Trotter factorization. This yields a hierarchy of approximations, one for each approximate factorization. For the purpose of a first benchmarking, we examine the approximate densities for a few cases with known exact densities and observe a very satisfactory, and encouraging, performance. As a bonus, we also obtain a simple fourth-order leapfrog algorithm for the symplectic integration of classical equations of motion.

Massively parallel symplectic algorithm for coupled magnetic spin dynamics and molecular dynamics

A parallel implementation of coupled spin–lattice dynamics in the LAMMPS molecular dynamics package is presented. The approach is very general, and can be applied to simple ferromagnets, magnetic alloys, or amorphous magnetic materials. Equations of motion for both spin only and coupled spin–lattice dynamics are first reviewed, including a detailed account of how magneto-mechanical potentials can be used to perform a proper coupling between spin and lattice degrees of freedom. A symplectic numerical integration algorithm is then presented which combines the Suzuki–Trotter decomposition for non-commuting variables and conserves the geometric properties of the equations of motion. The numerical accuracy of the serial implementation was assessed by verifying that it conserves the total energy and the norm of the total magnetization up to second order in the timestep size. Finally, a very general parallel algorithm is proposed that allows large spin–lattice systems to be efficiently simulated on large numbers of processors without degrading its mathematical accuracy. Its correctness as well as scaling efficiency were tested for realistic coupled spin–lattice systems, confirming that the new parallel algorithm is both accurate and efficient.

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.

A differential variational approach for handling fluid–solid interaction problems via smoothed particle hydrodynamics

The fluid–solid interaction (FSI) problem is customarily solved by starting with the fluid dynamics component. One chooses an established computational fluid dynamics (CFD) method and subsequently embeds the solid phase dynamics within the CFD solution leading to either a monolithic or a staggered/co-simulated solution. The approach discussed here takes the opposite tack. We start with a differential variational framework to handle the solid phase; i.e., the multi-body dynamics problem in the presence of contact, friction, and bilateral kinematic constraints. The dynamics of the fluid phase, which is captured via smoothed particle hydrodynamics (SPH), is subsequently embedded into this framework in which the incompressibility attribute of the flow is enforced via kinematic constraint equations that involve SPH particles. The resulting monolithic FSI solution methodology relies on a half-implicit symplectic time integration method that uses a matrix-free iterative approach to solve a cone constrained quadratic optimization problem at each time step. This problem yields the contact forces, friction forces, boundary condition Lagrange multipliers, fluid–solid coupling terms, and bilateral constraint Lagrange multipliers. The solution of the optimization problem represents the computationally taxing component of the method. Large integration time steps, tight enforcement of incompressibility, a unified approach for handling the fluid and solid phases, and linear scaling are listed as the attractive attributes of the proposed method. The numerical experiments reported include three validation studies (incompressibility, dam break, and sloshing), a scaling analysis, and a tracked vehicle fording simulation.

Lyapunov and reversibility errors for Hamiltonian flows

We discuss the stability of a Hamiltonian system by comparing the standard Lyapunov error (LE) with the forward error (FE) due to a small random perturbation. We introduce also the reversibility error (RE) where the evolution is computed forward up to time t and backwards to t=0 in presence of noise. This procedure has been investigated in the case of symplectic maps, but it turns out that the results are simpler in the case of a noisy flow, in the limit of zero noise amplitude. Indeed the stochastic processes defined by the displacement of the noisy orbit at time t for FE, or at time 0 for RE after the evolution up to time t, satisfy linear Langevin equations, are Gaussian processes, and the errors are just their root mean square deviations. All the errors are expressed in terms of the fundamental matrix L(t) of the tangent flow and can be evaluated numerically using a symplectic integrator. Letting eL(t) be the Lyapunov error and eR(t) be the reversibility error a very simple relation holds eR2(t)=∫0teL2(s)ds. The integral relation is quite natural since the local errors due to a random perturbations accumulate during the evolution whereas for the Lyapunov case the error is introduced only at time zero and propagated. The plot of errors for initial conditions in a Poincaré section reflects the phase portrait, whereas in the action plane it allows to single out the resonance strips. We have applied the method to a 3D Hamiltonian model H=H0(J)+λV(Θ), where analytic estimates can be obtained for the single resonances from perturbation theory. This allows to inspect the double resonance structure where the single resonance strips intersect. We have also considered the Hénon–Heiles Hamiltonian to show numerically the equivalence of the errors apart from a shift of 1/2 in the power law exponent in the case of regular orbits. The reversibility error method (REM), previously introduced as the error due to round off in the symplectic integration, appears to be comparable with RE also for the models considered here.

Long-term adaptive symplectic numerical integration of linear stochastic oscillators driven by additive white noise

In this paper, we present an adaptive variable step size numerical scheme for the integration of linear stochastic oscillator equations driven by additive Brownian white noise. We first show that traditional adaptive schemes based on local error estimation destroy the long-time behavior of the underlying method. As a remedy, we extend the idea presented in Hairer and Soderlind (SIAM J. Sci. Comput. 26(6), 1838–1851 2005) to the stochastic setting and show that using step density control mechanisms based on time regularization and local error tracking, we are able to obtain numerical schemes which preserve the important qualitative features of the solution process such as symmetry, time reversibility, symplecticity, linear growth rate of the second moment, and infinite oscillation. Numerical experiments confirm the theoretical findings of the paper.

Degenerate variational integrators for magnetic field line flow and guiding center trajectories

Symplectic integrators offer many benefits for numerically approximating solutions to Hamiltonian differential equations, including bounded energy error and the preservation of invariant sets. Two important Hamiltonian systems encountered in plasma physics—the flow of magnetic field lines and the guiding center motion of magnetized charged particles—resist symplectic integration by conventional means because the dynamics are most naturally formulated in non-canonical coordinates. New algorithms were recently developed using the variational integration formalism; however, those integrators were found to admit parasitic mode instabilities due to their multistep character. This work eliminates the multistep character, and therefore the parasitic mode instabilities via an adaptation of the variational integration formalism that we deem “degenerate variational integration.” Both the magnetic field line and guiding center Lagrangians are degenerate in the sense that the resultant Euler-Lagrange equations are systems of first-order ordinary differential equations. We show that retaining the same degree of degeneracy when constructing discrete Lagrangians yields one-step variational integrators preserving a non-canonical symplectic structure. Numerical examples demonstrate the benefits of the new algorithms, including superior stability relative to the existing variational integrators for these systems and superior qualitative behavior relative to non-conservative algorithms.

Field induced transient current in one-dimensional nanostructure

Field-induced transient current in one-dimensional nanostructures has been studied by a model of an electron confined in a 1D attractive Gaussian potential subjected both to electrodes at the terminals and to an ultrashort pulsed oscillatory electric field with the central frequency ω and the FWHM pulse width Γ. The time-propagation of the electron wave packet has been simulated by integrating the time-dependent Schrödinger equation directly relying on the second-order symplectic integrator method. The transient current has been calculated as the flux of the probability density of the escaping wave packet emitted from the downstream side of the confining potential. When a static bias-field E0 is suddenly applied, the resultant transient current shows an oscillatory decay behavior with time followed by a minimum structure before converging to a nearly constant value. The ω-dependence of the integrated transient current induced by the pulsed electric field has shown an asymmetric resonance line-shape for large Γ while it shows a fringe pattern on the spectral line profile for small Γ. These observations have been rationalized on the basis of the energy-level structure and lifetime of the quasibound states in the bias-field modified confining potential obtained by the complex-scaling Fourier grid Hamiltonian method.

Variational and symplectic integrators for satellite relative orbit propagation including drag

Orbit propagation algorithms for satellite relative motion relying on Runge–Kutta integrators are non-symplectic—a situation that leads to incorrect global behavior and degraded accuracy. Thus, attempts have been made to apply symplectic methods to integrate satellite relative motion. However, so far all these symplectic propagation schemes have not taken into account the effect of atmospheric drag. In this paper, drag-generalized symplectic and variational algorithms for satellite relative orbit propagation are developed in different reference frames, and numerical simulations with and without the effect of atmospheric drag are presented. It is also shown that high-order versions of the newly-developed variational and symplectic propagators are more accurate and are significantly faster than Runge–Kutta-based integrators, even in the presence of atmospheric drag.

Some remarks on the complex J-symmetric eigenproblem

The eigenproblem for complex J-symmetric matrices HC=[ACD−AT], A,C=CT, D=DT∈Cn×n is considered. A proof of the existence of a transformation to the complex J-symmetric Schur form proposed in [21] is given. The complex symplectic unitary QR decomposition and the complex symplectic SR decomposition are discussed. It is shown that a QR-like method based on the complex symplectic unitary QR decomposition is not feasible here. A complex symplectic SR algorithm is presented which can be implemented such that one step of the SR algorithm can be carried out in O(n) arithmetic operations. Based on this, a complex symplectic Lanczos method can be derived. Moreover, it is discussed how the 2n×2n complex J-symmetric matrix HC can be embedded in a 4n×4n real Hamiltonian matrix.

Explicit symplectic geometric algorithms for quaternion kinematical differential equation

Solving quaternion kinematical differential equations U+0028 QKDE U+0029 is one of the most significant problems in the automation, navigation, aerospace and aeronautics literatures. Most existing approaches for this problem neither preserve the norm of quaternions nor avoid errors accumulated in the sense of long term time. We present explicit symplectic geometric algorithms to deal with the quaternion kinematical differential equation by modelling its time-invariant and time-varying versions with Hamiltonian systems and adopting a three-step strategy. Firstly, a generalized Euler U+02BC s formula and Cayley-Euler formula are proved and used to construct symplectic single-step transition operators via the centered implicit Euler scheme for autonomous Hamiltonian system. Secondly, the symplecticity, orthogonality and invertibility of the symplectic transition operators are proved rigorously. Finally, the explicit symplectic geometric algorithm for the time-varying quaternion kinematical differential equation, i.e., a non-autonomous and non-linear Hamiltonian system essentially, is designed with the theorems proved. Our novel algorithms have simple structures, linear time complexity and constant space complexity of computation. The correctness and efficiencies of the proposed algorithms are verified and validated via numerical simulations.

Stability criteria for wide binary stars harboring Oort Clouds

Context. In recent years, several numerical studies have been done in the field of the stability limit. Although, many of them included the analysis of asteroids or planets, is not possible to find in the literature any work on how the presence of a binary star could affect other possible configurations in a three-body problem. In order to develop this subject we consider other structures like Oort Clouds in wide binary systems. Regarding the existence of Oort Clouds in extrasolar systems there are recent works that do not reject its possible existence.Aim. The aim of this work is to obtain the stability limit for Oort Cloud objects considering different masses of the secondary star and zero and non-zero inclinations of the particles. We improve our numerical treatment getting a mathematical fit that allows us to find the limit and compare our results with other previous works in the field.Methods. We use a symplectic integrator to integrate binary systems where the primary star is m1 = 1 M⊙ and the secondary, m2, takes 0.25 M⊙ and 0.66 M⊙ in two sets of simulations S1 and S2. The orbital parameters of the secondary star were varied in order to study different scenarios. We also used two different integration times (one shorter than the other) and included the presence of 1000 to 10 000 massless particles in circular orbits to form the Oort Cloud. The particles were disposed in four different inclination planes to investigate how the presence of the binary companion could affect the stability limit.Results. Using the Maximum Eccentricity Method, emax, together with the critical semimajor axis acrit we found that the emax criteria could reduce the integration times to find the limit. For those cases where the particles were in inclined orbits we found that there are particle groups that survive the integration time with a high eccentricity. These particle groups are found for our two sets of simulations, meaning that they are independent of the secondary mass. We also find for the co-planar case that the numerical value of the stability limit for retrograde orbits is higher than those found for prograde orbits. These results are in agreement with several published studies. Finally, the results obtained in this work allow us to build a numerical expression depending of the mass ratio, e2 and ip to find acrit, which can be compared with other recent works in the field.

New Constraints on Gliese 876—Exemplar of Mean-motion Resonance

Abstract Gliese 876 harbors one of the most dynamically rich and well-studied exoplanetary systems. The nearby M4V dwarf hosts four known planets, the outer three of which are trapped in a Laplace mean-motion resonance. A thorough characterization of the complex resonant perturbations exhibited by the orbiting planets, and the chaotic dynamics therein, is key to a complete picture of the system’s formation and evolutionary history. Here we present a reanalysis of the system using 6 yr of new radial velocity (RV) data from four instruments. These new data augment and more than double the size of the decades-long collection of existing velocity measurements. We provide updated estimates of the system parameters by employing a computationally efficient Wisdom–Holman N-body symplectic integrator, coupled with a Gaussian process (GP) regression model to account for correlated stellar noise. Experiments with synthetic RV data show that the dynamical characterization of the system can differ depending on whether a white-noise or correlated-noise model is adopted. Despite there being a region of stability for an additional planet in the resonant chain, we find no evidence for one. Our new parameter estimates place the system even deeper into resonance than previously thought and suggest that the system might be in a low-energy, quasi-regular double apsidal corotation resonance. This result and others will be used in a subsequent study on the primordial migration processes responsible for the formation of the resonant chain.

GENERIC Integrators: Structure Preserving Time Integration for Thermodynamic Systems

Abstract Thermodynamically admissible evolution equations for non-equilibrium systems are known to possess a distinct mathematical structure. Within the GENERIC (general equation for the non-equilibrium reversible–irreversible coupling) framework of non-equilibrium thermodynamics, which is based on continuous time evolution, we investigate the possibility of preserving all the structural elements in time-discretized equations. Our approach, which follows Moser’s [1] construction of symplectic integrators for Hamiltonian systems, is illustrated for the damped harmonic oscillator. Alternative approaches are sketched.

Structure-Preserving Model-Reduction of Dissipative Hamiltonian Systems

Reduced basis methods are popular for approximately solving large and complex systems of differential equations. However, conventional reduced basis methods do not generally preserve conservation laws and symmetries of the full order model. Here, we present an approach for reduced model construction, that preserves the symplectic symmetry of dissipative Hamiltonian systems. The method constructs a closed reduced Hamiltonian system by coupling the full model with a canonical heat bath. This allows the reduced system to be integrated with a symplectic integrator, resulting in a correct dissipation of energy, preservation of the total energy and, ultimately, in the stability of the solution. Accuracy and stability of the method are illustrated through the numerical simulation of the dissipative wave equation and a port-Hamiltonian model of an electric circuit.

Explicit symplectic algorithms based on generating functions for relativistic charged particle dynamics in time-dependent electromagnetic field

Relativistic dynamics of a charged particle in time-dependent electromagnetic fields has theoretical significance and a wide range of applications. The numerical simulation of relativistic dynamics is often multi-scale and requires accurate long-term numerical simulations. Therefore, explicit symplectic algorithms are much more preferable than non-symplectic methods and implicit symplectic algorithms. In this paper, we employ the proper time and express the Hamiltonian as the sum of exactly solvable terms and product-separable terms in space-time coordinates. Then, we give the explicit symplectic algorithms based on the generating functions of orders 2 and 3 for relativistic dynamics of a charged particle. The methodology is not new, which has been applied to non-relativistic dynamics of charged particles, but the algorithm for relativistic dynamics has much significance in practical simulations, such as the secular simulation of runaway electrons in tokamaks.