Symplectic Integrators Research Articles

Geometric exponential integrators

In this paper, we consider exponential integrators for semilinear Poisson systems. Two types of exponential integrators are constructed, one preserves the Poisson structure, and the other preserves energy. Numerical experiments for semilinear Possion systems obtained by semi-discretizing Hamiltonian PDEs are presented. These geometric exponential integrators exhibit better long time stability properties as compared to non-geometric integrators, and are computationally more efficient than traditional symplectic integrators and energy-preserving methods based on the discrete gradient method.

Symplectic integrators with adaptive time step applied to runaway electron dynamics

Studying the dynamics of runaway electrons has theoretical and practical significance. As the system is highly relativistic, multi-scale and nonlinear, accurate and efficient numerical methods with long-term stability are necessary. In this paper, we develop symplectic methods with adaptive time step and discuss the choice of step-size functions in accordance with the simulation problems for runaway electrons. In the implementation, in order to explore the practical impact of runaway electron dynamics, we use the electromagnetic field and some parameters often used in the study of plasma problems. Numerical results show that the new derived symplectic methods with adaptive time step exhibit good invariant-preserving property and superior stability over longtime integration. Moreover, with appropriate adaptive technique, the numerical efficiency in simulations is improved apparently. They are illustrated in the numerical experiments.

Simulating Maxwell–Schrödinger Equations by High-Order Symplectic FDTD Algorithm

A novel symplectic algorithm is proposed to solve the Maxwell-Schr\"odinger (M-S) system for investigating light-matter interaction. Using the fourth-order symplectic integration and fourth-order collocated differences, Maxwell-Schr\"odinger equations are discretized in temporal and spatial domains, respectively. The symplectic finite-difference time-domain (SFDTD) algorithm is developed for accurate and efficient study of coherent interaction between electromagnetic fields and artificial atoms. Particularly, the Dirichlet boundary condition is adopted for modeling the Rabi oscillation problems under the semi-classical framework. To implement the Dirichlet boundary condition, image theory is introduced, tailored to the high-order collocated differences. For validating the proposed SFDTD algorithm, three-dimensional numerical studies of the population inversion in the Rabi oscillation are presented. Numerical results show that the proposed high-order SFDTD(4,4) algorithm exhibits better numerical performance than the conventional FDTD(2,2) approach at the aspects of accuracy and efficiency for the long-term simulation. The proposed algorithm opens up a promising way towards a high-accurate energy-conservation modeling and simulation of complex dynamics in nanoscale light-matter interaction.

Motion and Interaction of Envelope Solitons in Schrödinger Equation Simulated by Symplectic Algorithm

The expression of Gaussian envelope soliton in Schrodinger equations are given and proved in this paper. According to the characteristics of the Gauss envelope soliton, further proposed that the interaction between Gaussian envelope solitons exists in Schrodinger equation. The symplectic algorithm for solving Schrodinger equation is proposed after analysis characteristics of Schrodinger equation. First, the Schrodinger equation is transformed into the standard Hamiltonian canonical equation by separating the real and imaginary parts of wave function. Secondly, the symplectic algorithm is implemented by using the Euler center difference method for the canonical equation. The conserved quantity of symplectic algorithm is given, and the stability of symplectic algorithm is proved. The numerical simulation experiment was carried out on Schrodinger equation in Gauss envelope soliton motion and multi solitons interaction. The experimental results show that the proposed method is correct and the symplectic algorithm is effective.

Symplectic integration of PDEs using Clebsch variables

Many PDEs (Burgers' equation, KdV, Camassa-Holm, Euler's fluid equations,...) can be formulated as infinite-dimensional Lie-Poisson systems. These are Hamiltonian systems on manifolds equipped with Poisson brackets. The Poisson structure is connected to conservation properties and other geometric features of solutions to the PDE and, therefore, of great interest for numerical integration. For the example of Burgers' equations and related PDEs we use Clebsch variables to lift the original system to a collective Hamiltonian system on a symplectic manifold whose structure is related to the original Lie-Poisson structure. On the collective Hamiltonian system a symplectic integrator can be applied. Our numerical examples show excellent conservation properties and indicate that the disadvantage of an increased phase-space dimension can be outweighed by the advantage of symplectic integration.

Comparative Study of Wellhole Surrounding Rock under Nonuniform Ground Stress

The stress analysis of the wellhole surrounding rock and the regular failure of the wellhole has always been a concern for the well builders. Firstly, the Hamilton canonical equations are obtained by using the Hamiltonian variational principle in the sector domain, and the zero eigensolution and nonzero eigensolutions of the homogeneous equation are solved. According to the Hamiltonian operator matrix with the orthogonal eigenfunction system, the special solution form of the nonhomogeneous boundary condition equation is obtained. Then, according to the principle of the same coefficient being equal, the relationship equation between the direction eigenvalue and the angle coefficient is obtained, from which the specific expression of the special solution of the equation can be determined. Furthermore, the analytical solution of the wellhole surrounding rock problem under nonuniform ground stress is obtained by using the linear elastic accumulative principle. Finally, a concrete example is given to compare the finite element method and the symplectic algorithm. The results are consistent, which ensures the accuracy and the reliability of the symplectic algorithm. The relationship between the circumferential stress distribution around the hole and the lateral pressure coefficient is further analyzed.

Computational efficiency of numerical integration methods for the tangent dynamics of many-body Hamiltonian systems in one and two spatial dimensions

We investigate the computational performance of various numerical methods for the integration of the equations of motion and the variational equations for some typical classical many-body models of condensed matter physics: the Fermi-Pasta-Ulam-Tsingou (FPUT) chain and the one- and two-dimensional disordered, discrete nonlinear Schrodinger equations (DDNLS). In our analysis we consider methods based on Taylor series expansion, Runge-Kutta discretization and symplectic transformations. The latter have the ability to exactly preserve the symplectic structure of Hamiltonian systems, which results in keeping bounded the error of the system's computed total energy. We perform extensive numerical simulations for several initial conditions of the studied models and compare the numerical efficiency of the used integrators by testing their ability to accurately reproduce characteristics of the systems' dynamics and quantify their chaoticity through the computation of the maximum Lyapunov exponent. We also report the expressions of the implemented symplectic schemes and provide the explicit forms of the used differential operators. Among the tested numerical schemes the symplectic integrators $ABA864$ and $SRKN^a_{14}$ exhibit the best performance, respectively for moderate and high accuracy levels in the case of the FPUT chain, while for the DDNLS models $s9\mathcal{ABC}6$ and $s11\mathcal{ABC}6$ (moderate accuracy), along with $s17\mathcal{ABC}8$ and $s19\mathcal{ABC}8$ (high accuracy) proved to be the most efficient schemes.

Structure‐preserving local optimal control of mechanical systems

SummaryWhile system dynamics are usually derived in continuous time, respective model‐based optimal control problems can only be solved numerically, ie, as discrete‐time approximations. Thus, the performance of control methods depends on the choice of numerical integration scheme. In this paper, we present a first‐order discretization of linear quadratic optimal control problems for mechanical systems that is structure preserving and hence preferable to standard methods. Our approach is based on symplectic integration schemes and thereby inherits structure from the original continuous‐time problem. Starting from a symplectic discretization of the system dynamics, modified discrete‐time Riccati equations are derived, which preserve the Hamiltonian structure of optimal control problems in addition to the mechanical structure of the control system. The method is extended to optimal tracking problems for nonlinear mechanical systems and evaluated in several numerical examples. Compared to standard discretization, it improves the approximation quality by orders of magnitude. This enables low‐bandwidth control and sensing in real‐time autonomous control applications.

Symplectic integrators: T + V revisited and round-off reduced

Symplectic integrators separate a problem into parts that can be solved in isolation, alternately advancing these sub-problems to approximate the evolution of the complete system. Problems with a single, dominant mass can use mixed-variable symplectic (MVS) integrators that separate the problem into Keplerian motion of satellites about the primary, and satellite-satellite interactions. Here, we examine T+V algorithms where the problem is separated into kinetic T and potential energy V terms. T+V integrators are typically less efficient than MVS algorithms. This difference is reduced by using different step sizes for primary-satellite and satellite-satellite interactions. The T+V method is improved further using 4th and 6th-order algorithms that include force gradients and symplectic correctors. We describe three 6th-order algorithms, containing 2 or 3 force evaluations per step, that are competitive with MVS in some cases. Round-off errors for T+V integrators can be reduced by several orders of magnitude, at almost no computational cost, using a simple modification that keeps track of accumulated changes in the coordinates and momenta. This makes T+V algorithms desirable for long-term, high-accuracy calculations.

Symmetric and symplectic exponential integrators for nonlinear Hamiltonian systems

This letter studies symmetric and symplectic exponential integrators when applied to numerically computing nonlinear Hamiltonian systems. We first establish the symmetry and symplecticity conditions of exponential integrators and then show that these conditions are extensions of the symmetry and symplecticity conditions of Runge–Kutta methods. Based on these conditions, some symmetric and symplectic exponential integrators up to order four are derived. Two numerical experiments are carried out and the results demonstrate the remarkable numerical behaviour of the new exponential integrators in comparison with some symmetric and symplectic Runge–Kutta methods in the literature.

Numerical discretization of completely integrable nonlinear Hamiltonian systems

Completely integrable nonlinear Hamiltonian systems offer attractive models for novel accelerator systems that are being developed to push the frontiers of high intensity and high power beam applications. However, all known numerical methods used to discretize and simulate these systems spoil the integrability property of the original, time-continuous, theoretically integrable system. Since the hallmark of integrability is the existence of sufficiently many functionally independent first integrals in involution, we investigated several methods' performance in terms of invariant preservation over longer timescales. We conclude that, while novel structure preserving integrators would be ideal, symplectic methods under mild restrictions offer adequate invariant preservation. We especially highlight the very good performance of the Stormer-Verlet method at low accuracy and a modified Lobatto method at high accuracy for the Integrable Optics Test Accelerator (IOTA). However, we show that symplectic integrators are limited by numerically induced resonances, which make interpretation of some simulation results challenging. These limitations and workarounds are discussed.

Masses of the Kepler-419 planets from transit timing variations analysis

We perform dynamical fits to the Transit Timing Variations (TTVs) of Kepler-419b. The TTVs from 17 Kepler quarters are obtained from Holczer et al (2016, ApJS 225,9). The dynamical fits are performed using the MultiNest Bayesian inference tool, coupled to an efficient symplectic N-body integrator. We find that the existing TTV data alone are able to uniquely constrain the planetary masses of Kepler-419b and c. Our estimates are in a good agreement with previous mass determinations that combined different techniques and observations, such as TTVs, radial velocity measurements, and the photoeccentric effect. As expected, however, our mass estimates have larger uncertainty. We study the global stability of the system within the parameters uncertainties to discard possible unstable solutions. We conclude that our method applied to the Holczer et al. data can provide reliable determinations of the planetary parameters. In the fore-coming work, we will use Holczer et al. (2016b) data to determine planetary parameters for a large set of Kepler systems.

Stability of planets in triple star systems

Context. Numerous theoretical studies of the stellar dynamics of triple systems have been carried out, but fewer purely empirical studies that have addressed planetary orbits within these systems. Most of these empirical studies have been for coplanar orbits and with a limited number of orbital parameters. Aims. Our objective is to provide a more generalized empirical mapping of the regions of planetary stability in triples by considering both prograde and retrograde motion of planets and the outer star; investigating highly inclined orbits of the outer star; extending the parameters used to all relevant orbital elements of the triple’s stars and expanding these elements and mass ratios to wider ranges that will accommodate recent and possibly future observational discoveries. Methods. Using N-body simulations, we integrated numerically the various four-body configurations over the parameter space, using a symplectic integrator designed specifically for the integration of hierarchical multiple stellar systems. The triples were then reduced to binaries and the integrations repeated to highlight the differences between these two types of system. Results. This established the regions of secular stability and resulted in 24 semi-empirical models describing the stability bounds for planets in each type of triple orbital configuration. The results were then compared with the observational extremes discovered to date to identify regions that may contain undiscovered planets.

Novel symplectic integrators for the Klein–Gordon equation with space- and time-dependent mass

We consider the numerical time-integration of the non-stationary Klein–Gordon equation with position- and time-dependent mass. A novel class of time-averaged symplectic splitting methods involving double commutators is analyzed and 4th- and 6th-order integrators are obtained. In contrast with standard splitting methods (that contain negative coefficients if the order is higher than two), additional commutators are incorporated into the schemes considered here. As a result, we can circumvent this order barrier and construct high order integrators with positive coefficients and a much reduced number of stages, thus improving considerably their efficiency. The performance of the new schemes is tested on several examples.

A high accurate hamiltonian nodal position finite element method for spatial cable structures undergoing long-term large overall motion

This paper addresses the challenges faced by the error accumulation over long-term numerical calculation for dynamic modeling of spatial flexible cable structures undergoing large translational and rotational motion. A high accurate Hamiltonian nodal position finite element method is proposed to deal with the challenges. The new nodal position finite element discrete formulation is derived by Hamiltonian theory and Green strain theory with full expression of global stiffness matrices without additional simplifications. Symplectic difference algorithm is built for numerical solution to optimally preserve the energy, momenta and area (volume) of the phase space. Forth-order closed Newton-Cotes numerical integration is applied to calculate the aerodynamic drag force. The Symplectic conservation feature of the proposed method is validated by the dynamics of a long-period classical pendulum. The numerical accuracy and stability of the proposed method are validated by LS-DYNA simulations for a flexible polyethylene rubber conical pendulum, the experiments of a three-dimensional circularly towed cable and the experiments of a free swing cable. The present algorithm is compared with the conventional second-order Runge-Kutta algorithm. All validations and comparisons indicate that the proposed method is stable and highly accurate.

Estimation of Surface Fault Displacement by High Performance Computing

Reliable estimation of surface fault displacements is crucial to the safety of nuclear power plant facilities. It is necessary to develop a numerical method for the estimation. In the study, we develop a finite element method in which the following two functions are implemented: (1) a symplectic time integration of an explicit scheme to properly conserve the energy of the system; and (2) rigorously formulated joint elements of high order. The finite element method is enhanced with parallel computing capability. We apply the developed method to solve simple three-dimensional models of faults embedded in a rock mass. It includes a comparison of results from quasi-static and dynamic simulations and investigation of the sensitivity of results to the shear stiffness on faults. In the study, we propose capacity computing with a quasi-static simulation for uncertainty quantification.

Symplectic orbit propagation based on Deprit’s radial intermediary

The constantly challenging requirements for orbit prediction have opened the need for better onboard propagation tools. Runge-Kutta (RK) integrators have been widely used for this purpose; however RK integrators are not symplectic, which means that RK integrators may lead to incorrect global behavior and degraded accuracy. Emanating from Deprit’s radial intermediary, obtained by the elimination of the parallax transformation, we present the development of symplectic integrators of different orders for spacecraft orbit propagation. Through a set of numerical simulations, it is shown that these integrators are more accurate and substantially faster than Runge-Kutta-based methods. Moreover, it is also shown that the proposed integrators are more accurate than analytic propagation algorithms based on Deprit’s radial intermediary solution, and even other previously-developed symplectic integrators.

Symplectic integration of boundary value problems

Symplectic integrators can be excellent for Hamiltonian initial value problems. Reasons for this include their preservation of invariant sets like tori, good energy behaviour, nonexistence of attractors, and good behaviour of statistical properties. These all refer to {\em long-time} behaviour. They are directly connected to the dynamical behaviour of symplectic maps $\varphi\colon M\to M$ on the phase space under iteration. Boundary value problems, in contrast, are posed for fixed (and often quite short) times. Symplecticity manifests as a symplectic map $\varphi\colon M\to M'$ which is not iterated. Is there any point, therefore, for a symplectic integrator to be used on a Hamiltonian boundary value problem? In this paper we announce results that symplectic integrators preserve bifurcations of Hamiltonian boundary value problems and that nonsymplectic integrators do not.

Dynamical cartography of Earth satellite orbits

We have carried out a numerical investigation of the coupled gravitational and non-gravitational perturbations acting on Earth satellite orbits in an extensive grid, covering the whole circumterrestrial space, using an appropriately modified version of the SWIFT symplectic integrator, which is suitable for long-term (∼120 years) integrations of the non-averaged equations of motion. Hence, we characterize the long-term dynamics and the phase-space structure of the Earth-orbiter environment, starting from low altitudes (∼400 km) and going up to the GEO region and beyond. This investigation was done in the framework of the EC-funded “ReDSHIFT” project, with the purpose of enabling the definition of passive debris removal strategies, based on the use of physical mechanisms inherent in the complex dynamics of the problem (i.e., resonances). Accordingly, the complicated interactions among resonances, generated by different perturbing forces (i.e., lunisolar gravity, solar radiation pressure, tesseral harmonics in the geopotential) are accurately depicted in our results, where we can identify the regions of phase space where the motion is regular and long-term stable and regions for which eccentricity growth and even instability due to chaotic behavior can emerge. The results are presented in an “atlas” of dynamical stability maps for different orbital zones, with a particular focus on the (drag-free) range of semimajor axes, where the perturbing effects of the Earth’s oblateness and lunisolar gravity are of comparable order. In some regions, the overlapping of the predominant lunisolar secular and semi-secular resonances furnish a number of interesting disposal hatches at moderate to low eccentricity orbits. All computations were repeated for an increased area-to-mass ratio, simulating the case of a satellite equipped with an on-board, area-augmenting device. We find that this would generally promote the deorbiting process, particularly at the transition region between LEO and MEO. Although direct reentry from very low eccentricities is very unlikely in most cases of interest, we find that a modest “delta-v” (ΔV) budget would be enough for satellites to be steered into a relatively short-lived resonance and achieve reentry into the Earth’s atmosphere within reasonable timescales (∼50 years).

GpuSPHASE—A shared memory caching implementation for 2D SPH using CUDA (new version announcement)

Smoothed particle hydrodynamics (SPH) is a meshless Lagrangian method that has been successfully applied to computational fluid dynamics (CFD), solid mechanics and many other multi-physics problems. gpuSPHASE is a graphics processing unit (GPU) accelerated solver for 2D SPH simulations that is optimized for relatively small numbers of particles but long physical simulation times as required to solve transport phenomena in process engineering applications. The software aims at a low latency execution pipeline with thousands of iterations per second on contemporary GPU accelerator cards, while using symplectic time integration for long term stability of the simulations. For this, a novel caching algorithm for CUDA shared memory is proposed and implemented. gpuSPHASE is validated against SPHERIC test cases and the performance is evaluated and compared to state of the art implementations of the SPH method. New version program summaryProgram Title: gpuSPHASEProgram Files doi:http://dx.doi.org/10.17632/8vpwwh8th5.2Licensing provisions: GPLv3Programming language: CUDA C++Journal reference of previous version: Comput. Phys. Comm., 213 (2017), 165–180Does the new version supersede the previous version?: YesReasons for the new version: The software is in active development and thus several new features have been added since the initial release. Those improvements are related to accuracy, maintainability and platform independence.Summary of revisions:•Build process: The previous distribution was relatively complicated to build due to the external dependencies of the software. This process has been optimized such that CMake is used for all platforms to build the software. For performance reasons several decisions have to be taken at compile time, which can be configured through the CMake configuration. External dependencies are platform independently retrieved and included using CMake features.•Code: While optimized for performance, zero-cost decorators were added to improve the readability of the code. Furthermore, a weak form of type punning was used previously, which was entirely replaced by a custom defined type to eliminate the possibility of undefined behavior.•Floating point precision: As described in the literature and in the original version of the manuscript, 32 bit floating point precision is generally sufficient for SPH simulations of water [1,2]. While not optimized for the opportunistic caching strategy, 64 bit floating point precision has been added. Due to the significant performance loss of consumer grade GPUs for FP64 computations the performance may be considerably worse. FP64 precision can be used to check if FP32 is sufficient and thus detect simulations where FP64 is required. This feature does not influence the already implemented cell relative position calculations that allow for large computational domains without any loss in accuracy [2]. Improvements can be observed e.g. in the Poiseuille flow test case released with DualSPHysics 4.2 [3] with the configuration of ν=10−6 [m2.s−1], ρ=1000 [kg.m−3] and a body force of Fx=10−4 [m.s−2], as shown in Fig. 1.•Particle shifting: The diffusion based particle shifting method proposed by Skillen et al. (2013) [4] is implemented. As suggested by Lind et al. (2012) [5], the particle concentration is augmented with the factor fij to prevent pairing instability for free surface flows [6]. The improvements are best observed with Poiseuille flow configuration of ν=10−2 [m2.s−1], ρ=1 [kg.m−3] and a body force of Fx=10−1 [m.s−2] [7]. This test case does not converge with gpuSPHASE using single or double precision due to the problems described by [8]. Enabling particle shifting allows to execute the simulation until t>100, as shown in Fig. 2.•Advective-diffusion and heat transfer: Heat transfer is modeled based on the temperature T and the thermal conductivity k as dT∕dt=ρ−1cp−1k∇2T. The SPH discretization is implemented according to Alshaer et al. (2017) [9] and successfully verified [10]. The mass transfer of concentration C is modeled using the advective-diffusion equation dC∕dt=D∇2C−∇⋅(vC). The model is implemented according to the discretization of Aristodemo et al. [11] and successfully verified [12].Nature of problem: gpuSPHASE aims at free surface fluid dynamics simulations of long running physical phenomena that must be calculated in the order of real-time. Due to the long physical time and thus the requirement of millions to billions of iterations, the code is restricted to 2D. This restriction allows for optimizations that would not be possible for more generic approaches [13].Solution method: gpuSPHASE is a 2D SPH solver for CUDA capable devices that is optimized for the computation of real-time simulations. Compile time definitions allow to optimize the code for fast execution of the CUDA kernel functions. Symplectic time integration and correction methods enable stable simulations of long physical time problems.