Symplectic Numerical Integration Research Articles

Strong backward error analysis of symplectic integrators for stochastic Hamiltonian systems

Backward error analysis is a powerful tool in order to detect the long-term conservative behavior of numerical methods. In this work, we present a long-term analysis of symplectic stochastic numerical integrators, applied to Hamiltonian systems with multiplicative noise. We first compute and analyze the associated stochastic modified differential equations. Then, suitable bounds for the coefficients of such equations are provided towards the computation of long-term estimates for the Hamiltonian deviations occurring along the aforementioned numerical dynamics. This result generalizes Benettin-Giorgilli Theorem to the scenario of stochastic symplectic methods. Finally, specific numerical methods are considered, in order to provide a numerical evidence confirming the effectiveness of the theoretical investigation.

Symplectic numerical methods in optics and imaging: ray tracing in spherical gradient-index lenses and computer-generated image rendering.

This paper provides an introduction to symplectic numerical integration techniques and examines various optical applications. We first outline the fundamentals of Hamiltonian optics and detail the construction of a symplectic method via the splitting technique. Numerical experiments involving a selection of spherically symmetric gradient-index lenses compare the accuracy of various first-, second-, and fourth-order symplectic methods with equivalent nonsymplectic methods. The best-performing methods are then further tested as part of an image rendering task involving nonlinear ray tracing, comparing the trace time required by each method. Future improvements, recommendations, and uses for symplectic ray tracing are also considered.

A new symplectic integrator for stochastic Hamiltonian systems on manifolds

This paper presents a significant development in the field of geometric numerical stochastic integration schemes. Specifically, the Geometric Symplectic Itô-Taylor 1.0 strong numerical integration scheme, tailored for Hamiltonian systems evolving on 2-sphere or S2 manifold is proposed. The core of this proposed algorithm centers around a new weak symplectic condition. This condition ensures the numerical stability for long duration time integration that can be readily achieved through numerical steps rather than relying solely on analytical satisfaction of the condition. This proposed advancement over the state-of-the-art caters to improved accuracy and stability of numerical simulations in manifold-based dynamical systems for long duration simulations. Through extensive analysis and numerical experiments, the effectiveness and reliability of the proposed scheme and symplectic criterion are validated. The findings of this study promises to offer valuable insights for researchers and practitioners working in the field of symplectic numerical integration methods.

Symplectic numerical integration for Hamiltonian stochastic differential equations with multiplicative Lévy noise in the sense of Marcus

In this paper, we propose a symplectic numerical integration method for a class of Hamiltonian stochastic differential equations with multiplicative Lévy noise in the sense of Marcus. We first construct a general symplectic Euler scheme for these equations, then we prove its convergence. In addition, we provide realizable numerical implementations for the proposed symplectic Euler scheme in detail. Some numerical experiments are conducted to demonstrate the effectiveness and superiority of the proposed method by the simulations of its orbits, Hamiltonian and convergence order over a long time interval. The results show the applicability of the methods considered.

Structure-preserving techniques in accelerator physics

To very good approximation, particularly for hadron machines, charged-particle trajectories in accelerators obey Hamiltonian mechanics. During routine storage times of eight hours or more, such particles execute some revolutions about the machine, oscillations about the design orbit, and passages through various bending and focusing elements. Prior to building, or modifying, such a machine, we seek to identify accurately the long-term behaviour and stability of particle orbits over such large numbers of interactions. This demanding computational effort does not yield easily to traditional methods of symplectic numerical integration, including both explicit Yoshida-type and implicit Runge–Kutta or Gaussian methods. As an alternative, one may compute an approximate one-turn map and then iterate that map. We describe some of the essential considerations and techniques for constructing such maps to high order and for realistic magnetic field models. Particular attention is given to preserving the symplectic condition characteristic of Hamiltonian mechanics.

Symplectic integration with non-canonical quadrature for guiding-center orbits in magnetic confinement devices

We study symplectic numerical integration of mechanical systems with a Hamiltonian specified in non-canonical coordinates and its application to guiding-center motion of charged plasma particles in magnetic confinement devices. The technique combines time-stepping in canonical coordinates with quadrature in non-canonical coordinates and is applicable in systems where a global transformation to canonical coordinates is known explicitly but its inverse is not. A fully implicit class of symplectic Runge-Kutta schemes has recently been introduced and applied to guiding-center motion by Zhang et al. (2014) [9]. Here a generalization of this approach with emphasis on semi-implicit partitioned schemes is described together with methods to enhance performance, in particular avoiding evaluation of non-canonical variables at full time steps. For application in toroidal plasma confinement configurations with nested magnetic flux surfaces a global canonicalization of coordinates for the guiding-center Lagrangian by a spatial transform is presented that allows for pre-computation of the required map in a parallel algorithm in the case of time-independent magnetic field geometry. Guiding-center orbits are studied in stationary magnetic equilibrium fields of an axisymmetric tokamak and a realistic three-dimensional stellarator configuration. Superior long-term properties of symplectic methods are demonstrated in comparison to a conventional adaptive Runge-Kutta scheme. Finally statistics of fast fusion alpha particle losses over their slowing-down time are computed in the stellarator field on a representative sample, reaching a speed-up of the symplectic Euler scheme by more than a factor three compared to usual Runge-Kutta schemes while keeping the same statistical accuracy and linear scaling with the number of computing threads.

Discrete-time port-Hamiltonian systems: A definition based on symplectic integration

We introduce a new definition of discrete-time port-Hamiltonian (PH) systems, which results from structure-preserving discretization of explicit PH systems in time. We discretize the underlying continuous-time Dirac structure with the collocation method and add discrete-time dynamics by the use of symplectic numerical integration schemes. The conservation of a structural discrete-time energy balance – expressed in terms of the discrete-time Dirac structure – extends the notion of symplecticity of geometric integration schemes to open systems. We discuss the energy approximation errors in the context of the presented definition and show that their order for linear PH systems is consistent with the order of the numerical integration scheme. Implicit Gauss–Legendre methods and Lobatto IIIA/IIIB pairs for partitioned systems are examples for integration schemes that are covered by our definition. The statements on the numerical energy errors are illustrated by elementary numerical experiments.

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.

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.

Molecular dynamics with rigid bodies: Alternative formulation and assessment of its limitations when employed to simulate liquid water.

Sets of atoms collectively behaving as rigid bodies are often used in molecular dynamics to model entire molecules or parts thereof. This is a coarse-graining strategy that eliminates degrees of freedom and supposedly admits larger time steps without abandoning the atomistic character of a model. In this paper, we rely on a particular factorization of the rotation matrix to simplify the mechanical formulation of systems containing rigid bodies. We then propose a new derivation for the exact solution of torque-free rotations, which are employed as part of a symplectic numerical integration scheme for rigid-body dynamics. We also review methods for calculating pressure in systems of rigid bodies with pairwise-additive potentials and periodic boundary conditions. Finally, simulations of liquid phases, with special focus on water, are employed to analyze the numerical aspects of the proposed methodology. Our results show that energy drift is avoided for time step sizes up to 5 fs, but only if a proper smoothing is applied to the interatomic potentials. Despite this, the effects of discretization errors are relevant, even for smaller time steps. These errors induce, for instance, a systematic failure of the expected equipartition of kinetic energy between translational and rotational degrees of freedom.

On the application of improved symplectic integrators in Hamiltonian Monte Carlo

ABSTRACTWe explore the construction of new symplectic numerical integration schemes to be used in Hamiltonian Monte Carlo and study their efficiency. Integration schemes from Blanes et al., and a new scheme are considered as candidates to the commonly used leapfrog method. All integration schemes are tested within the framework of the No-U-Turn sampler (NUTS), both for a logistic regression model and a student t-model. The results show that the leapfrog method is inferior to all the new methods both in terms of asymptotic expected acceptance probability for a model problem and the efficient sample size per computing time for the realistic models.

Hamiltonian B-series and a Lie algebra of non-rooted trees

The conditions for a B-series to be a Hamiltonian vector field imply that it may be given as a series indexed by free trees, i.e. trees without root. At present, the pre-Lie structure of rooted trees plays an important role in the study of numerical methods for ordinary differential equations, as does the associated Lie bracket on rooted trees obtained by antisymmetrization. We give a substitute for this Lie bracket defined on free trees that reflects the Lie bracket of Hamiltonian B-series, and illustrate an application of this to the backward error analysis of symplectic numerical integrators.

Bill2d — A software package for classical two-dimensional Hamiltonian systems

We present Bill2d, a modern and efficient C++ package for classical simulations of two-dimensional Hamiltonian systems. Bill2d can be used for various billiard and diffusion problems with one or more charged particles with interactions, different external potentials, an external magnetic field, periodic and open boundaries, etc. The software package can also calculate many key quantities in complex systems such as Poincaré sections, survival probabilities, and diffusion coefficients. While aiming at a large class of applicable systems, the code also strives for ease-of-use, efficiency, and modularity for the implementation of additional features. The package comes along with a user guide, a developer’s manual, and a documentation of the application program interface (API). Program summaryProgram title: Bill2dCatalogue identifier: AEYL_v1_0Program summary URL:http://cpc.cs.qub.ac.uk/summaries/AEYL_v1_0.htmlProgram obtainable from: CPC Program Library, Queen’s University, Belfast, N. IrelandLicensing provisions: GNU General Public License, version 3No. of lines in distributed program, including test data, etc.: 37098No. of bytes in distributed program, including test data, etc.: 1155037Distribution format: tar.gzProgramming language: C++(14).Computer: Tested on x86 and x86 64 architectures.Operating systems: Tested on Linux, and OS X versions 10.9–10.11.Has the code been vectorized or parallelized?: Shared memory parallelization when simulating ensembles of systems. Vectorization of operations with R2 vectors.RAM: Simulation dependent: kilobytes to gigabytesClassification: 4.3, 7.8, 7.9, 7.10, 16.9.External routines: Boost, CMake, GSL, HDF5; and optionally Google-Mock, GoogleTest, and DoxygenNature of problem: Numerical propagation of classical two-dimensional single and many-body systems, possibly in a magnetic field, and calculation of relevant quantities such as Poincaré sections, survival probabilities, diffusion co-efficients, etc.Solution method: Symplectic numerical integration of Hamilton’s equations of motion in Cartesian coordinates, or solution of Newton’s equations of motion if in a magnetic field. The program implements several well-established algorithms.Restrictions: Pointlike particles with equal masses and charges, although the latter restrictions are easy to lift.Unusual features: Program is efficient, extremely modular and easy to extend, and allows arbitrary particle–particle interactions.Additional comments: The source code is also available at https://bitbucket.org/solanpaa/bill2d. See README for locations of user guide, developer manual, and API docs.Running time: From milliseconds to days, depends on type of simulation.

Heisenberg-Integrable Spin Systems

We investigate certain classes of integrable classical or quantum spin systems. The first class is characterized by the recursively defined property P saying that the spin system consists of a single spin or can be decomposed into two uniformly coupled or disjoint subsystems with property P. For these systems the time evolution can be explicitly calculated. The second class consists of spin systems where all non-zero coupling constants have the same strength (spin graphs) possessing N − 1 independent, commuting constants of motion of Heisenberg type. These systems are shown to have the above property P and can be characterized as spin graphs not containing chains of length four as vertex-induced sub-graphs. We completely enumerate and characterize all spin graphs up to N = 5 spins. Applications to the construction of symplectic numerical integrators for non-integrable spin systems are briefly discussed.

Dynamics and stability of telluric planets within the habitable zone of extrasolar planetary systems

Aims. We study gravitational perturbation effects of observed giant extrasolar planets on hypothetical Earth-like planets in the context of the three-body problem. This paper considers a large parameter survey of different orbital configuration of two extrasolar giant planets (HD 70642b and HD 4208b) and compares their dynamical effect on Earth-mass planetary orbits initially located within the respective habitable terrestrial region. We are interested in determining giant-planet orbit (and mass) parameters that favor the condition to render an Earth-mass planet to remain on a stable and bounded orbit within the continuous habitable zone.Methods. We applied symplectic numerical integration techniques to studying the short and long term time evolution of hypothetical Earth-mass planets that are treated as particles. In addition, we adopt the MEGNO technique to obtain a complete dynamical picture of the terrestrial phase space environment. Both multi-particle and single-particle simulations were performed to follow an Earth-mass planet in the habitable region and its subsequent long term evolution.Results. Our numerical simulations show that giant planets should be on circular orbits to minimize the perturbative effect on terrestrial orbits. The orbit eccentricity (and hence proximity) is the most important orbital parameter of dynamical significance. The most promising candidate for maintaining an Earth-mass planet on a stable and bounded orbit well-confined to the continuous habitable zone is HD 70642b. Even the large planetary mass of HD 70642b renders an Earth-mass planet habitable during the complete lifetime of the host star. The results allow us to extrapolate similar observed systems and points the necessity further constraining the uncertainty range in giant planet orbital eccentricity by future follow-up observations.

High-Precision, Symplectic Numerical, Relative Orbit Propagation

This paper presents a numerical method to propagate relative orbits. It can handle an arbitrary number of zonal and tesseral terms in the geopotential. This method relies on defining a relative Hamiltonian, which describes both the absolute and the relative motion of two satellites. The solution is separated into an analytical Keplerian part and a symplectic numerical integration part. The algorithm is designed to conserve the constants of the motion, resulting in better long-term accuracy. We present results for a broad range of scenarios with large separations and show that submeter accuracy is possible over five days of propagation, with a geopotential model containing 36 terms in tesseral and zonal harmonics. These results are valid for eccentricities reaching 0.5. Furthermore, the relative propagation scheme is significantly faster than differencing two absolute orbit propagations.

Symplectic tracing of high‐energy charged particles in the inner magnetosphere

To investigate the radiation belt, the Gauss Runge‐Kutta method (RKG), one of the symplectic numerical integration methods, is introduced and tested by tracing protons and electrons in geomagnetic dipole field. The tracings of protons having the energies from 10 keV to 10 MeV with pitch angles from 90° to 30° at L ∼ 5 are practiced by fourth‐ and sixth‐order RKG with time step 0.016 of the Larmor period on the equatorial plane. The accuracy of calculations are examined with Cartesian, cylindrical, and polar coordinates. With all coordinates, we confirmed that the calculations are stable for 10,000,000 steps with RKG. On the other hand, the calculations for protons in the same conditions with standard Runge‐Kutta methods fail within 10,000 steps. The tracings of relativistic electrons having energies from 10 to 150 MeV with pitch angles from 90° to 30° at L ∼ 5 are also practiced. For the electrons above 50 MeV, the calculations are successful for 10,000,000 steps with the time step 0.1 × RE/c by fourth and sixth RKG with cylindrical coordinates. The calculations by standard Runge‐Kutta method are failed in any coordinate system, and the calculations by fourth or sixth RKG with Cartesian or polar coordinates are also failed for all energy range above. The calculations for the electrons below 10 MeV with this time step size are all failed by any methods and with any coordinate systems mentioned above. For the effective calculations, for example, for 10 MeV electron with pitch angle 90° by fourth order RKG, we must use shorter timescale size than 0.2 of the above size.

The stability of the orbits of Earth-mass planets in the habitable zone of 47 Ursae Majoris

We have investigated whether Earth-mass planets could survive in the habitable zone (HZ) of the 47 Ursae Majoris system. Mixed-variable symplectic numerical integration has been used to investigate the orbits of putative Earth-mass planets. Whereas the 47 UMa system as previously known, with just one giant planet, could have Earth-mass planets that remain confined to the HZ for a fairly wide range of initial obital parameters, the second (outer) giant now known to be present has reduced the range significantly. There are however confined orbits particularly if either the present eccentricity of the outer giant is close to zero, or the giant masses are close to the minimum values. We have also shown that the eccentricity of the outer giant's orbit is unlikely to exceed about 0.15 at minimum giant mass and 0.13 at 1.5 times this minimum.

The resonance theory for stochastic layers in nonlinear dynamic systems

A theory of stochastic layers is developed for a better understanding of the resonant mechanism of stochastic layers in nonlinear Hamiltonian systems. A criterion based on an accurate whisker map and resonant conditions is developed for prediction of the onset of resonance in the stochastic layer. The onset of a specific primary resonance between the periodic forcing and periodic orbit of the integrable Hamiltonian system in the stochastic layer is predicted analytically. A forced, twin-well Duffing oscillator is investigated as a sample problem for prediction of a specified resonance in the stochastic layer. Verification of the analytical prediction is carried out through a symplectic numerical integration scheme. The analytical and numerical results are in good agreement.

Multi-symplectic integrators: numerical schemes for Hamiltonian PDEs that conserve symplecticity

The symplectic numerical integration of finite-dimensional Hamiltonian systems is a well established subject and has led to a deeper understanding of existing methods as well as to the development of new very efficient and accurate schemes, e.g., for rigid body, constrained, and molecular dynamics. The numerical integration of infinite-dimensional Hamiltonian systems or Hamiltonian PDEs is much less explored. In this Letter, we suggest a new theoretical framework for generalizing symplectic numerical integrators for ODEs to Hamiltonian PDEs in R 2 : time plus one space dimension. The central idea is that symplecticity for Hamiltonian PDEs is directional: the symplectic structure of the PDE is decomposed into distinct components representing space and time independently. In this setting PDE integrators can be constructed by concatenating uni-directional ODE symplectic integrators. This suggests a natural definition of multi-symplectic integrator as a discretization that conserves a discrete version of the conservation of symplecticity for Hamiltonian PDEs. We show that this approach leads to a general framework for geometric numerical schemes for Hamiltonian PDEs, which have remarkable energy and momentum conservation properties. Generalizations, including development of higher-order methods, application to the Euler equations in fluid mechanics, application to perturbed systems, and extension to more than one space dimension are also discussed.