Symplectic Integrators Research Articles

Step-size effect in the time-transformed leapfrog integrator on elliptic and hyperbolic orbits

ABSTRACT A drift-kick-drift (DKD) type leapfrog symplectic integrator applied for a time-transformed separable Hamiltonian (or time-transformed symplectic integrator; TSI) has been known to conserve the Kepler orbit exactly. We find that for an elliptic orbit, such feature appears for an arbitrary step size. But it is not the case for a hyperbolic orbit: When the half step size is larger than the conjugate momenta of the mean anomaly, a phase transition happens and the new position jumps to the non-physical counterpart of the hyperbolic trajectory. Once it happens, the energy conservation is broken. Instead, the kinetic energy minus the potential energy becomes a new conserved quantity. We provide a mathematical explanation for such phenomenon. Our result provides a deeper understanding of the TSI method, and a useful constraint of the step size when the TSI method is used to solve the hyperbolic encounters. This is particular important when an (Bulirsch–Stoer) extrapolation integrator is used together, which requires the convergence of integration errors.

Spin-dependent transient current in transistor-like nanostructures

Transient current in transistor-like nanostructures has been studied by a model of a few electrons confined in a one-dimensional effective potential consisting of three quantum wells, ‘source’, ‘gate’, and ‘drain’. The time-dependent Schrödinger equation for the electrons has been integrated relying on the symplectic integrator method and the transient current has been calculated as the flux of the probability density of electrons absorbed by the complex absorbing potential placed at the far edge of the drain region. The electrons are initially placed in the source domain as their lowest-energy state for a given spin multiplicity and the source–drain current has been calculated for different gate potential heights. The current for different spin configurations has shown strong emission at different values of the gate potential, suggesting use of the studied nanostructures for extracting current with a specific spin configuration from spin-unpolarized normal current. Dependence of the current emission on electron correlation has also been studied by changing the size of the source domain. The current has shown appreciable differences for different spin configurations for the medium and strong confinement regimes, while these differences become smaller for smaller confinement and tend to diminish in the weak limit of confinement. This observed trend has been rationalized on the basis of the formation of the Wigner lattice states.

The use of Kepler solver in numerical integrations of quasi-Keplerian orbits

ABSTRACT A Kepler solver is an analytical method used to solve a two-body problem. In this paper, we propose a new correction method by slightly modifying the Kepler solver. The only change to the analytical solutions is that the obtainment of the eccentric anomaly relies on the true anomaly that is associated with a unit radial vector calculated by an integrator. This scheme rigorously conserves all integrals and orbital elements except the mean longitude. However, the Kepler energy, angular momentum vector, and Laplace–Runge–Lenz vector for perturbed Kepler problems are slowly varying quantities. However, their integral invariant relations give the quantities high-precision values that directly govern five slowly varying orbital elements. These elements combined with the eccentric anomaly determine the desired numerical solutions. The newly proposed method can considerably reduce various errors for a post-Newtonian two-body problem compared with an uncorrected integrator, making it suitable for a dissipative two-body problem. Spurious secular changes of some elements or quasi-integrals in the outer Solar system may be caused by short integration times of the fourth-order Runge–Kutta algorithm. However, they can be eliminated in a long integration time of 108 yr by the proposed method, similar to Wisdom–Holman second-order symplectic integrator. The proposed method has an advantage over the symplectic algorithm in the accuracy but gives a larger slope to the phase error growth.

Random and interval uncertain symplectic methods for linear Birkhoff equations and their comparison study

With the increasing complexity and design requirement of engineering structures, the numerical computation plays an increasingly important role in the structural analysis. The symplectic algorithm for Birkhoffian systems has advantages over the traditional algorithm in accuracy and stability, and general mechanical systems can be transformed into Birkhoffian systems in theory. Additionally, the influence of uncertainties in practical engineering cannot be ignored. However, the research foundation of symplectic algorithms for Birkhoffian systems considering uncertainties is still weak. In this paper, the random and interval uncertain symplectic algorithms for linear Birkhoffian systems are proposed, respectively, which can be applied to the numerical calculation of the structural dynamic response with uncertainties. Furthermore, the compatibility between the two methods is studied. Based on the small parameter perturbation theory, the uncertain perturbation equation of the uncertain linear Birkhoffian system is established. According to different types of uncertainties, based on the probability theory and interval mathematical method, expressions of the mean and variance of random Birkhoffian systems, as well as the median value and radius of interval Birkhoffian systems are derived. By combining aforementioned methods with Birkhoffian symplectic algorithms, the random and interval perturbation methods for calculating the uncertain response of linear Birkhoffian systems are proposed, respectively. Since the uncertainty exists objectively, the uncertain features of the response should not be incompatible when characterization and analysis methods for uncertainties are different. Based on the Chebyshev inequality, this paper deduces the probabilistic and interval boundaries from the result calculated by proposed methods and the comparison between the two kinds of boundaries. It is mathematically proved that the probabilistic upper bound is less than the interval upper bound and the probabilistic lower bound is larger than the interval lower bound. Finally, examples are given to demonstrate the effectiveness of the proposed method in calculating the response of uncertain dynamics and its superiority over the traditional numerical methods.

Mass-zero constrained molecular dynamics for electrode charges in simulations of electrochemical systems.

Classical molecular dynamics simulations have recently become a standard tool for the study of electrochemical systems. State-of-the-art approaches represent the electrodes as perfect conductors, modeling their responses to the charge distribution of electrolytes via the so-called fluctuating charge model. These fluctuating charges are additional degrees of freedom that, in a Born-Oppenheimer spirit, adapt instantaneously to changes in the environment to keep each electrode at a constant potential. Here, we show that this model can be treated in the framework of constrained molecular dynamics, leading to a symplectic and time-reversible algorithm for the evolution of all the degrees of freedom of the system. The computational cost and the accuracy of the new method are similar to current alternative implementations of the model. The advantage lies in the accuracy and long term stability guaranteed by the formal properties of the algorithm and in the possibility to systematically introduce additional kinematic conditions of arbitrary number and form. We illustrate the performance of the constrained dynamics approach by enforcing the electroneutrality of the electrodes in a simple capacitor consisting of two graphite electrodes separated by a slab of liquid water.

Symplectic Computations of Fast Ion Trajectory and Radial Diffusion Coefficient in MHD Perturbed Tokamak

The magneto-hydrodynamics (MHD) modes existing in conducting plasmas can greatly affect charged particle’s dynamics. Low frequency MHD perturbations, called the neo-classical tearing (NTM) modes, give rise to magnetic islands which can deteriorate fast ion confinement, especially in tokamak plasmas. Here, a Hamiltonian guiding-center code, based on symplectic algorithm which allows efficient computations of energetic ion trajectories and transport coefficients, has been applied to magnetically perturbed tokamak plasmas. The effects of NTM perturbations on the energetic ion trajectories and diffusion coefficient have been presented. It is demonstrated that the inclusion of such steady state magnetic perturbations in the guiding center drift orbit computations escalates the loss of energetic ions via radial diffusion. Moreover, the broader NTMs lead to a significant enhancement of the radial diffusion coefficient.

A symplectic integrator for molecular dynamics on a hypersphere

We present a reversible and symplectic algorithm called ROLL, for integrating the equations of motion in molecular dynamics simulations of simple fluids on a hypersphere $\mathcal{S}^d$ of arbitrary dimension $d$. It is derived in the framework of geometric algebra and shown to be mathematically equivalent to algorithm RATTLE. An application to molecular dynamics simulation of the one component plasma is briefly discussed.

Lie Group Cohomology and (Multi)Symplectic Integrators: New Geometric Tools for Lie Group Machine Learning Based on Souriau Geometric Statistical Mechanics.

In this paper, we describe and exploit a geometric framework for Gibbs probability densities and the associated concepts in statistical mechanics, which unifies several earlier works on the subject, including Souriau’s symplectic model of statistical mechanics, its polysymplectic extension, Koszul model, and approaches developed in quantum information geometry. We emphasize the role of equivariance with respect to Lie group actions and the role of several concepts from geometric mechanics, such as momentum maps, Casimir functions, coadjoint orbits, and Lie-Poisson brackets with cocycles, as unifying structures appearing in various applications of this framework to information geometry and machine learning. For instance, we discuss the expression of the Fisher metric in presence of equivariance and we exploit the property of the entropy of the Souriau model as a Casimir function to apply a geometric model for energy preserving entropy production. We illustrate this framework with several examples including multivariate Gaussian probability densities, and the Bogoliubov-Kubo-Mori metric as a quantum version of the Fisher metric for quantum information on coadjoint orbits. We exploit this geometric setting and Lie group equivariance to present symplectic and multisymplectic variational Lie group integration schemes for some of the equations associated with Souriau symplectic and polysymplectic models, such as the Lie-Poisson equation with cocycle.

Optimized Lie–Trotter–Suzuki decompositions for two and three non-commuting terms

Lie–Trotter–Suzuki decompositions are an efficient way to approximate operator exponentials exp(tH) when H is a sum of n (non-commuting) terms which, individually, can be exponentiated easily. They are employed in time-evolution algorithms for tensor network states, digital quantum simulation protocols, path integral methods like quantum Monte Carlo, and splitting methods for symplectic integrators in classical Hamiltonian systems. We provide optimized decompositions up to order t6. The leading error term is expanded in nested commutators (Hall bases) and we minimize the 1-norm of the coefficients. For n=2 terms, several of the optima we find are close to those in McLachlan (1995). Generally, our results substantially improve over unoptimized decompositions by Forest, Ruth, Yoshida, and Suzuki. We explain why these decompositions are sufficient to efficiently simulate any one- or two-dimensional lattice model with finite-range interactions. This follows by solving a partitioning problem for the interaction graph.

GFiRe—Gauge Field integrator for Reheating

We present a new numerical algorithm and code, GFiRe, for solving the non-linear evolution of Abelian gauge fields coupled to complex scalar fields in homogeneous and isotropic spacetimes. We adopt a hybrid approach to solving the system: the spatial derivatives are discretized using standard Lattice Gauge Field Theory techniques, whereas the time evolution of the fields and scale factor is implemented with explicit, composite, symplectic integrators. An important property of our compound algorithm is that the discretized Gauss constraint is respected exactly, regardless of the order of the symplectic integrator. This remains true even when the background expansion is computed “self-consistently”; that is, when the expansion history is computed using spatial averaged components of the energy momentum tensor in the simulation volume. Hence, our code can also be used in cases where the fields dominate the energy density of the universe, for example, during reheating after inflation. We test the algorithm in scenarios of reheating where the inflaton is a complex scalar field with a potential ∝(2|φ|2−v2)2 and is coupled to an Abelian gauge field. Tracing the evolution of the system through complex dynamics (including resonant excitation of fields, backreaction, formation of solitons, and changes in the equation of state) in a self-consistently expanding universe, we find the energy conservation violation (<10−4) to be very stable and the Gauss constraint violation (<10−6) to be dominated by differencing noise.

Accelerated methods for direct computation of fusion alpha particle losses within, stellarator optimization

Accelerated statistical computation of collisionless fusion alpha particle losses in stellarator configurations is presented based on direct guiding-centre orbit tracing. The approach relies on the combination of recently developed symplectic integrators in canonicalized magnetic flux coordinates and early classification into regular and chaotic orbit types. Only chaotic orbits have to be traced up to the end, as their behaviour is unpredictable. An implementation of this technique is provided in the code SIMPLE (symplectic integration methods for particle loss estimation, Albert et al., 2020b, doi:10.5281/zenodo.3666820). Reliable results were obtained for an ensemble of 1000 orbits in a quasi-isodynamic, a quasi-helical and a quasi-axisymmetric configuration. Overall, a computational speed up of approximately one order of magnitude is achieved compared to direct integration via adaptive Runge–Kutta methods. This reduces run times to the range of typical magnetic equilibrium computations and makes direct alpha particle loss computation adequate for use within a stellarator optimization loop.

Preservation of Bifurcations of Hamiltonian Boundary Value Problems Under Discretisation

We show that symplectic integrators preserve bifurcations of Hamiltonian boundary value problems and that nonsymplectic integrators do not. We provide a universal description of the breaking of umbilic bifurcations by nonysmplectic integrators. We discover extra structure induced from certain types of boundary value problems, including classical Dirichlet problems, that is useful to locate bifurcations. Geodesics connecting two points are an example of a Hamiltonian boundary value problem, and we introduce the jet-RATTLE method, a symplectic integrator that easily computes geodesics and their bifurcations. Finally, we study the periodic pitchfork bifurcation, a codimension-1 bifurcation arising in integrable Hamiltonian systems. It is not preserved by either symplectic on nonsymplectic integrators, but in some circumstances symplecticity greatly reduces the error.

Optimal Third-Order Symplectic Integration Modeling of Seismic Acoustic Wave Propagation

ABSTRACTSeismic wavefield modeling based on the wave equation is widely used in understanding and predicting the dynamic and kinematic characteristics of seismic wave propagation through media. This article presents an optimal numerical solution for the seismic acoustic wave equation in a Hamiltonian system based on the third-order symplectic integrator method. The least absolute truncation error analysis method is used to determine the optimal coefficients. The analysis of the third-order symplectic integrator shows that the proposed scheme exhibits high stability and minimal truncation error. To illustrate the accuracy of the algorithm, we compare the numerical solutions generated by the proposed method with the theoretical analysis solution for 2D and 3D seismic wave propagation tests. The results show that the proposed method reduced the phase error to the eighth-order magnitude accuracy relative to the exact solution. These simulations also demonstrated that the proposed third-order symplectic method can minimize numerical dispersion and preserve the waveforms during the simulation. In addition, comparing different central frequencies of the source and grid spaces (90, 60, and 20 m) for simulation of seismic wave propagation in 2D and 3D models using symplectic and nearly analytic discretization methods, we deduce that the suitable grid spaces are roughly equivalent to between one-fourth and one-fifth of the wavelength, which can provide a good compromise between accuracy and computational cost.

An efficient numerical treatment for the asymptotic behaviour of the nonlinear Airy-type problems

This study focuses on symplectic integrators for numerical evaluation of the asymptotic solutions of the nonlinear Airy-type equations obtained by reducing the nonlinear dispersive equations. Since the nature of Airy-type equations has both highly oscillatory slow decay and exponential fast decay, most of classical integrators are not able to correctly exhibit challenging physical behaviour. We use specially designed symplectic integrators combining splitting methods with Magnus integrators to catch asymptotic behaviour of nonlinear Airy-type equations efficiently, even for large step sizes. Efficiency of the proposed methods for given problems is discussed. Moreover, numerical results obtained by the proposed methods are compared with the existing results in the literature.

A slow-down time-transformed symplectic integrator for solving the few-body problem

ABSTRACTAn accurate and efficient method dealing with the few-body dynamics is important for simulating collisional N-body systems like star clusters and to follow the formation and evolution of compact binaries. We describe such a method which combines the time-transformed explicit symplectic integrator and the slow-down method. The former conserves the Hamiltonian and the angular momentum for a long-term evolution, while the latter significantly reduces the computational cost for a weakly perturbed binary. In this work, the Hamilton equations of this algorithm are analysed in detail. We mathematically and numerically show that it can correctly reproduce the secular evolution like the orbit averaged method and also well conserve the angular momentum. For a weakly perturbed binary, the method is possible to provide a few orders of magnitude faster performance than the classical algorithm. A publicly available code written in the c++ language, sdar, is available on github. It can be used either as a standalone tool or a library to be plugged in other N-body codes. The high precision of the floating point to 62 digits is also supported.

Symplectic ray tracing based on Hamiltonian optics in gradient-index media.

A method based on Hamiltonian optics for ray tracing through gradient-index (GRIN) media is proposed. The ray equation that describes light-ray paths can be written in the form of the Hamiltonian equations. Although the Hamiltonian equations can be numerically calculated using a finite-difference explicit method, deviations from the exact equations are generally inevitable at subsequent time steps. An optical Hamiltonian can be constructed of two independent terms, i.e., one term dependent on position and the other term dependent on momentum. The symplectic integrator is applicable to such a separable optical Hamiltonian system and makes the optical Hamiltonian equations form invariant at each time step of numerical calculations. Accuracies of light-ray paths calculated using the first-order symplectic ray tracing in GRIN lenses approximate those calculated on the basis of the fourth-order Runge-Kutta algorithm, which shows the promising potential of the symplectic-ray-tracing method.

Application of PDS–FEM to simulate dynamic crack propagation and supershear rupture

Dynamic crack propagation problems, including intersonic rupture, are simulated with two Hamiltonian based formulations of particle discretization scheme (PDS) FEM: the traditional displacement momentum form and the strain momentum form, for which consistent momentum conserving and symplectic integration schemes are derived. Numerical results are verified, and validated by comparing with photoelastic observations of a dynamic mode-I crack captured with a 1 Mfps camera. Both methods are successful in accurately reproducing the crack patterns observed in classical 2D and 3D dynamic rupture scenarios, as well as the near crack tip stress field during the propagation. The two methods are found to be numerically indifferentiable, although the displacement based method offers a significantly better computational performance. As a demonstrative application, we simulate the super-shear rupture in earthquakes, modeling the contact at the fault surface by a linear slip weakening friction law. The Burridge–Andrews mechanism naturally appears in the simulations, making the crack front jump from the sub-Rayleigh regime to the intersonic regime and propagate while producing shear Mach cones.

An 𝒪(n) framework for internal coordinate molecular dynamics applicable to molecules with arbitrary constraints and geometries

ABSTRACTFor molecules with constraints such as fixed lengths and angles, it is more efficient to consider the molecular movement in the space of generalised internal coordinate than in Cartesian coordinate. This paper presents a new framework in the simulation of molecular movement, especially for macro-molecules with massive length and angle constraints. The generalised forces are calculated to invert the dense mass matrix for integrating the constraints equation. The inverting of the mass matrix was made based on distance descending ordering method. The method does an reordering of the internal variables to make the Cholesky decomposition need no fill-in, which promises the time complexity in doing mass matrix inverting. The method was extended for application to loop structures. It is found that the mass matrix would become singular when the bond angles approach 0 or π; thus a rotation convention switch method was proposed to resolve the singularity. The time complexity has been demonstrated and the length and angle constraints can be arbitrarily applied. The long-time energy conservation in NVE ensemble was compared for results from symplectic and non-symplectic time integrators. The non-symplectic fourth-order Runge–Kutta method still has satisfactory long-time energy conservation if using small time step.

Discrete-Time Control Design Based on Symplectic Integration: Linear Systems

We propose a novel approach to design discrete-time state feedback controllers for sampled control systems with guaranteed stability under arbitrary sampling times that fulfill the Nyquist-Shannon condition. The key idea is borrowed from symplectic integration and backward error analysis of Hamiltonian systems: The closed-loop target system is the discretized version of a continuous-time system with appropriately shaped Hamiltonian, where a symplectic scheme is used for discretization. We adopt this argumentation for a systematic discrete-time design procedure for controllable linear systems based on the implicit midpoint rule. We motivate the approach on the example of two basic linear systems under zero order hold sampling, we show the construction of target systems with desired eigenvalues based on the discrete-time controller canonical form, and we illustrate the quality of the main result on a random sixth order system.

Forward Analysis of GPR for Underground Pipes Using CUDA-Implemented Conformal Symplectic Euler Algorithm

Ground-penetrating radar (GPR) is widely used in the detection and positioning of underground facilities. Through the inversion analysis of echo signals of GPR, information such as pipe material, burial depth and location of pipelines can be obtained. Unfortunately, underground pipelines are cylindrical, the ladder approximation method used in traditional forward models produces certain errors. In this study, an accurate and efficient numerical model of GPR forward model in underground pipelines is established using symplectic Euler algorithm, graphics processing unit (GPU) acceleration technology and surface conformal technology. With a Ricker wavelet pulse as the GPR source, the convolution perfectly matched layer (CPML) is incorporated in the symplectic Euler algorithm and shown to be effective to truncate the symplectic Euler computational domain. Through the simulation study of different underground pipeline models, GPR image characteristics of the metal pipeline, plastic pipeline and concrete pipeline filled with air and water are obtained. According to the numerical simulation results, parallel conformal symplectic Euler algorithm effectively reduces the false diffracted waves caused by ladder approximation and improves the computational efficiency of the model in metallic and non-metallic media.