Introduction

The Hamiltonian system plays a vital part in describing the evolution of a physical system. Various numerical schemes, including the Euler, Runge–Kutta, linear multistep methods, and some low-order Taylor methods, have been designed to describe the dynamical system. However, these methods are non-structure-preserving. In addition, they lead to unstable computation or incorrect solutions. The symplectic method [1–3], a common structure-preserving method, conserves the area or volume of the system during computation. The square conservation also preserves the structure through conserving the length of the simulated system [4]. These structure-preserving methods allow the Hamiltonian constant to remain constant, or only vary periodically during the entire integration time range [0,t]. The structure-preserving methods have the advantages in that they provide a true long-term trajectory of the simulated system and stabilize the computation process.

However, these structure-preserving methods need to be improved. First, when tackling the nonlinear Hamiltonian systems, some of the symplectic methods strongly depend on the generating function. These implicit methods are applied to solve nonlinear algebraic equations at each step, and thus efficiency becomes a problem. Some symplectic methods are based on the Runge–Kutta method, the order of which is generally less than 10 (rarely more than 15), to avoid a complicated procedure. Some other high-order explicit symplectic methods are used to study separable Hamiltonian systems [2,5], but these explicit methods usually have an order of less than 10, and are limited to separable Hamiltonian systems.

Second, although symplectic methods have the advantage of choosing a larger step size than in a classical approach and thus saving simulation time, such an increase in step-size causes large errors in the primary variables. This occurs despite that there is no modification in the trajectory structure. Gladman [6] has said that “…the conservation of the integrals is not a problem for the symplectic integration algorithms (SIAs) but the phase errors can still be uncomfortable after a large number of orbits. If one wanted to integrate the Earth for the lifetime of the solar system it is doubtful that these two SIAs could perform the ~10⁹ orbit integration reliably. This is not necessarily too disheartening, since no other integration scheme known to the authors could perform the integration either.”

The phase trajectory of a Hamiltonian system is one of the most basic requirements. However, given the period of the Hamiltonian system, the symplectic method provides no special insight, and only gives approximate numerical periods with the precision proportional to the order of method and step-size. The long-term integration of a dynamical system is a challenge but urgent task in many fields. Orders of magnitude of time periods in physics range from the Planck time (5.39 × 10⁻⁴⁴ s) to the age of universe (about 10¹⁷ s). Thus, a meaningful non-dimensional time for a dynamical system is within 10⁶⁰ orders of magnitude. Simulation of the position of dynamical variables (not the trajectory structure) in the ultra long-term (i.e., t = 10⁶⁰) is still a time consuming task, even for a periodic dynamical system. The purpose of this study is to provide a reliable scheme to solve the issues mentioned above.

Materials and Methods

The parallel multiple-precision Taylor (PMT) method [7–12] is originally designed to solve nonlinear chaotic systems. It provides ultra-high reliable solutions for longer times than other methods. The order of the PMT method can be very high compared to other traditional approaches. Here, the application of the PMT method to a nonlinear Hamiltonian system is examined. The orbits of the system for two atoms with Morse potential energy [13] are (1) The initial values are x₀ = 0 and , where x is the displacement and p is the momentum of the particles. Using the substitution y = e^−x yields (2) where y₀ = 1 and , and is a constant.

The Taylor series expansions relevant to solving (2) are (3) where the coefficients are given by and ; and h is the step-size. The coefficients are determined from the initial conditions , and the relations (4)

More details associated with the parallel scheme are referred to Wang et al. [7]. The solution of (1) may be expressed as x_n+1 = −ln y_n+1, after obtaining the values of y_n+1 and p_n+1. If the order, M, is larger than 100, the parallel scheme greatly reduces the computation time; if M is smaller than 30, one CPU is generally sufficient to perform the calculation on a reasonable timescale.

The symplectic method used to solve (1) is a 2^nd order explicit method (SE2), as discussed by Qin et al. [14]: where c₁ = 0, c₂ = 1, , f(x) = −(e^−2x − e^−x), and g(p) = p.

The second and more complex example is a non-separable Hamiltonian system, which is often broadly defined as: (5) The Hamiltonian is , and the initial values are q₀ = 0, p₀ = 1: where b(t) = sin q, g(t) = cos q, c(t) = p sin q, d(t) = p + cos q and the k^th Taylor coefficients are b_k, g_k, c_k, and d_k. Therefore, other coefficients are: , , , , , d_m = α_m + g_m, and the initial coefficients are α₀ = p₀, β₀ = q₀, b₀ = sin q₀, g₀ = cos q₀, c₀ = p₀ sin q₀, and d₀ = p₀ + cos q₀.

The conservation of the Hamiltonian indicates . In the present numerical simulation, if |ΔH| ≤ 10⁻¹⁶ (i.e., the smallest relative error in double-precision floating point arithmetic); it is regarded as unchanged. When H is a non-zero constant, is a criterion for a large Hamiltonian.

For instance, there is a Hamiltonian constant, H = −0.01, for Eq 2 when the numerical solutions of p and x are p_N and x_N, respectively. The error of the Hamiltonian is

Since and guarantees |ΔH| < 10⁻¹⁶. The numerical error at time t indicates that |p_N − p| ≤ C₁h^M+1 and |y_N − y| ≤ C₂h^M+1. where C₁ and C₂ are constants, and . Since |p| and |y| are bounded variables, where C is a constant that satisfies |p|C₁ + |y|C₂ + C₂ ≤ C.

With a step-size of h = 0.01, a high enough order M is chosen to guarantee . In practice, the order of M can be easily determined by several numerical runs without knowing the value of C. By using this high-order method, the structure-preserving solution of the original equation is obtained by numerical means. In fact, because it is very easy to increase M with the Taylor series method, we can make |ΔH| even smaller to meet Eq 2.

In the direct simulation of Eq 2 with t = 10⁷, a 20-order PMT scheme is used to achieve the simulation results. Fig 1 illustrations of the direct simulation of Eq 2. In Fig 1C, the PMT method is shown to predict the correct trajectory structure (x-p plane) along with a correct cycle of x (Fig 1A). During the entire computation time range, the Hamiltonian H approaches a constant (Fig 1D), while Fig 1E shows that H varies periodically and has larger errors when using the SE2 method. From Fig 1B we could found a more important issue is that the error in x increases as the simulation time increases. Thus, the position of x is not reliable at times longer than 10⁵ time units.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 1.

Illustrations of the direct simulation of Eq 2: (a) variable x by PMT; (b) error of x computed by the SE2 method; (c) structure of the x-p plane by PMT; (d) the Hamiltonian H by PMT to T = 10⁷; (e) error of the Hamiltonian H by the SE2 method. Panels (a–c) have step-size h = 0.01; (d–e) have step-size h = 0.1.

https://doi.org/10.1371/journal.pone.0163303.g001

Fig 2 is the direct simulation of Eq 5 by the 20-order PMT method. Meanwhile, using the PMT method to solve the non-separable equation, the variable (Fig 2B) and the Hamiltonian H (Fig 2A) are simulated well. The Hamiltonian H remains to be a constant throughout the simulation in Fig 2A.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 2.

The direct simulation of Eq 5 by the 20-order PMT method, with a step-size of h = 0.01: (a) Hamiltonian H, (b) variable p versus time (the first 20 time units).

https://doi.org/10.1371/journal.pone.0163303.g002

These results are all indicates that the PMT scheme can be used to simulate dynamical system with very high precision. Furthermore, the PMT scheme is a self-verifying scheme, as discussed in [7]. This verification scheme is a standard operation for PMT and CNS [12,15,16] numerical experiments.

The (Forward period analysis) FPA computation procedure is combined by two stages, one is the period finding and another is application of the founded period to do long term simulation. The FPA (and PMT scheme) are all suitable to solve the linear problem, but here we use a nonlinear system as objection to test the scheme’s performance. The more details of FPA are reported in the following section.

Results

FPA stage 1: The period discover of a Hamiltonian system

Establishing the phase trajectories of Hamiltonian systems is a basic requirement, which can be achieved by the symplectic method as well as PMT. However, the symplectic method only gives approximate numerical periods, with a precision proportional to the order of the method and step-size.

The key of numerical methods to identify the period of a dynamical system is to find out when the solutions return to the initial values (if the system is defined by n-th 1^rst order differential equation, the n variables must return to their initial values simultaneously). The integration time between the start point and the repeat point is thus approximately the period. Generally, the period obtained by numerical method in this way varies. The error between a variable (such as x) and its corresponding repeat point (defined as E_x, for example E_x ≤ 10⁻³⁰) indicates the uncertainty (Δt) in estimating the period. Because E_x is small and , Δt ⋅ p ≈ E_x. The standard division of p is obtained through numerical experiments, , as an averaged value of p such that Δt ∼ E_x / σ(p). This formula suggests that the error bounds of a typical period have a magnitude of 10⁻³⁰.

The forward period analysis (FPA) method is proposed to obtain the period for Eq 2. The first stage of FPA is a pre-computation to find a suitable residual interval. The computation starts with y = 1, and , and a time-step size of h₀ = 0.01. At each step, the values of y_k and y_k−1 are checked to find out the approximately first period. The first repeat position satisfies y ≈ 1 and , determined within the interval of [T_l,T_h], where T_l is the lower and T_h the higher bound. This interval is the time it takes for y to cross the base line y = 1 (from the y > 1 to the y < 1 direction, the steps that reach the lower bound are defined as k). The first period is now between T_l = kh₀ and T_h = (k + 1)h₀ time units, and thus the period T ≈ T_l. The error of the period is about (h₀ + Δt) ≈ 10⁻², and thus the precision of this forward period analysis method is mostly dependent on the last computation step-size h_F (in this stage, h_F = h₀).

The second stage is the post-computation at the residual interval [kh₀,(k + 1)h₀] where k = Int(t / T) is a positive integer number. This interval can be separated into subintervals by the dichotomy method. Denoting the whole interval as before, [T_l,T_h], with T_l = kh₀ and T_h = (k + 1)h₀, a new step size h_F = h_F / 2 is chosen to separate the interval into [T_l,T_l + h_F] and [T_l + h_F,T_h]. If the value of y at T_l + h_F does not cross the base line, then T_l = T_l + h_F; otherwise, T_h = T_l + h_F and then the operation is repeated in the new interval [T_l,T_h]. The dichotomy method maintains h_F smaller than the magnitude of E_T / σ(p), and thus total error of the period is dominated by ΔT (ΔT < T_h − T_l ≈ 10⁻³⁰). The computation cost for the dichotomy method in the last interval is about 30log₂ 10 ≈ 100 loops, while that in the pre-computation stage is about T / 0.01 loops. The value of T can be roughly estimated from Fig 3A (T is within a 45 time unit). The total loops for obtaining the period of Eq 2 are within 5000 loops.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 3.

Demonstration of the FPA method to obtain the period for Eq 2: (a) the numerical result of variable y by the PMT method in the interval [0,100]; and (b) variable y in the enlarged interval of [43,45].

https://doi.org/10.1371/journal.pone.0163303.g003

The above procedures are also suitable for SE2 and other symplectic methods; if we choose the step-size (h) for SE2, the period with precision at the magnitude of h² is obtained. Applying the dichotomy method at the last interval for SE2, the error at T_l and T_h must first be confirmed to be small enough, and this is not the superiority for SE2. Since a decreas in h greatly increases the computation time, if a more accurate period is required, for example ΔT = 10⁻³⁰, h ≈ 10⁻¹⁵ must be set, and this requires 10¹⁵ loops to finish the computation.

In addition, applying self-verification requires reducing step-size to be about h/100 or less in SE2 to guarantee that the reference solution is more accurate than the solution in which step-size is h, and this requires 100 times more loops than the computation process. While self-verification of the PMT method only requires increasing the order M to a bigger value for example (M+10), this does not increase the number of computation loops of the verification process, but only the time cost per loop. The increasing time cost in one integration loop is insignificant when M is large (for example M>100). As a consequence, the PMT method is efficiently verified.

Fig 3 is the demonstration of the FPA method to obtain the period for Eq 2. After we obtain the roughly value of period T from Fig 3A. we then can obtain a more precise value from the enlarged time axis in Fig 3B. In Fig 3B the first stage of FPA is to determine the residual interval of Eq 2 as [44.42,44.43], i.e. T_l = 44.42 and T_h = 44.43. The FPA procedure in this interval, and the values of 30 and 100 significant digits, are listed in Table 1. To the author’s knowledge, this level of accuracy has never been reported.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 1. The period of a Morse system obtained by FPA with M = 200, h = 0.01, and the precision we use is 2000 bits.

https://doi.org/10.1371/journal.pone.0163303.t001

FPA stage 2: The application of FPA in long-term simulations

Before demonstrating an ultra long-term simulation of Eq 2, a simple periodic dynamical system is analyzed to determine the most important parameters in the long-term computation. The simple dynamical system is defined by (6) and the initial values are

The issue is how to obtain 16 significant digits of x(t) at t = 10³⁰. Since the analytical solution of this equation is the result should be x(10³⁰) = sin(10³⁰). As sin(t − 2πk) = sin(t), the result is given by sin(10³⁰) = sin(10³⁰ − 2πk), and , i.e., the integer part of . The precision of sin(10³⁰) is dependent on the precision of the approximation to 2π. The reference value of π with 50 significant digits and the computed k are listed in Table 2, and the double-precision (16 significant digits) results are also compared. From Table 2, note that the k values corresponding to the two different precisions of π are different. The different k-values cause the residual of t, i.e., T_r = 10³⁰ − 2πk, to be more uncertain. Therefore, precision to 16 significant digits for sin(10³⁰) is not possible in a double-precision float platform.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 2. The values of π and k.

https://doi.org/10.1371/journal.pone.0163303.t002

The above example indicates that the reliable long-term computation of a periodic system is dependent on the precision of the period; 2π can be regarded as the period in this example. The relative error bound, ε, is estimated for the period to limit the computation error at t to Δx = 10⁻¹⁶. The true residual time can be written as and the residual time induced by numerical error is

The first restriction of ε guarantees that and are the same. Under this situation, we have T_r = 10³⁰ − 2πk, and the difference between T_r and is Δt = 2πkε. The error bounds, ε, satisfying 2πkε < 10⁻¹⁶ will guarantee Δx < 10⁻¹⁶. Since 0 < 10³⁰ − 2πk < 2π, 2πk ≈ 10³⁰ and ε < 10⁻¹⁶/2πk = 10⁻⁴⁶ is the relative error bound of 2π. The 50 significant digits of π satisfy this error bound. Thus, T_r = 3.231831977487846 and sin(10³⁰) = −0.090116901912138058.

The analysis of error bounds for a general periodic dynamical system is the same as the example by changing the period from 2π to T. Thus, ε < E(t)/t is the period error bound, where t is the simulation time and E(t) the required relative error bound for the output. The fast computation of sin(10³⁰) is a benefit from the pre-known of precise value of π. However, for a general periodic dynamical system such as Eq 2 there is no such pre-knowledge for the period T, hence FPA is proposed to obtain the precise T first.

Long-term simulation by FPA is achieved by dividing the long-term (t ≫ T) computation into two parts: one is the detection of period of a cycle; the other is the simulation of the residual time, equal to t − kT. Because the computation error of the symplectic method generally propagates linearly with time [17,18], such that an increase of t times requires a t^1/M times smaller step-size to control the error, where M denotes the order of the symplectic method. The computation complexity for time t in unit loops is O(t^1+1/M). However, applying the FPA procedure with PMT helps to reduce the computation from O(t^1+1/M) to O(T / h₀) in the first stage, and to O(ln t) in the second stage, i.e., O(ln t + T / h₀).

As illustrated in Table 3, the long-term computation of a dynamical system uses at least O(t) loops without FPA, but with FPA this is cut to O(ln t + T / h₀) at most. This improvement greatly decreases the CPU time cost and makes many unsolvable long-term problems be reliably solvable.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 3. The loops used in the different schemes.

https://doi.org/10.1371/journal.pone.0163303.t003

Given the period obtained by FPA (see Table 1), the variable values for the Morse system will be computed quickly. The results of selected times from (10 to 10⁶⁰) are listed in Table 4, with the corresponding values obtained by direct integration. The results of direct integrations with FPA and PMT agree well. The results of the SE2 method have errors from the early integration stage, and the results beyond 10⁵ are incorrect. It took about one day to obtain the result at t = 10⁷ by direct integration, so obtaining a result at 10⁶⁰ is a seemingly impossible task for direct integration, but with the help of FPA, reliable results can be obtained.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 4. The variable p obtained by FPA and direct simulation for a Morse system.

https://doi.org/10.1371/journal.pone.0163303.t004

Another classical dynamical system is the mathematical pendulum. The Hamiltonian of a pendulum system is , and the dynamical equation is (7) with initial values

The period of this system approaches to 2π, while the initial momentum p → 0 [17]. Table 5 lists the period corresponding to the initial condition, q = 0, with different momenta, p. All periods are accurate to 50 significant digits. As illustrated in Table 5, the period approaches 2π when p decreases from 1 to 10⁻³⁰. This experiment again proves the correctness of FPA.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 5. The period of a mathematical pendulum system obtained by FPA with M = 200, h = 0.01, and the precision we use is 2000 bits.

https://doi.org/10.1371/journal.pone.0163303.t005

It is very fast applying FPA to obtain a period for these demonstration systems, and the computation can be finished within 1 minute on a Linux system with an Intel Xeon 2.5 Ghz CPU. The long-time scope solutions can be obtained within 1 minute by FPA.

Discussion and Conclusion

Using the FPA, we obtain the periods of some classical Hamiltonian systems, with the accuracy of 100 significant digits. We also confirm that reliable solutions within the time range t ∈ [0,10⁶⁰] can be obtained. To the best of the author’s knowledge, this accuracy of time period for long-term solutions of Hamiltonian systems has not been reported before. The FPA method provides a powerful tool to gain time-effective ultra long-term reliable solutions of periodic systems.

The FPA procedure works well in conjunction with the traditional symplectic method and the PMT method. Generally, symplectic methods with different orders require different subroutines to conduct the computation. In contrast, it is relatively easy to change the order of the Taylor series method, so that it provides the flexibility to carry out simulations with different orders of accuracy for one system. The PMT method reasonably simulates the Morse system for t ∈ [0,10⁷], but for much longer times simulation it is hard without FPA. This is demonstrated using a simple example. For more complex systems, higher order PMT approaches can be used. Indeed, details of an example for application of a high-order PMT method in direct simulationg of the Henon–Heiles system is referred to Liao [15].

The Taylor series method has a good convergence property when the order is high enough [9]. This feature can enlarge the step-size h to 0.01 for Eq 2, and increase the simulation speed. The result obtained by the Taylor series method not only maintains ΔH ≃ 0, but also reduces numerical errors. The PMT method is not a structure-preserving method, but it can preserve the structure well by numerical means. Consequently, it can be used as an alternative of symplectic methods for the computation of simple Hamiltonian systems. Moreover, PMT can be applied to some non-separable nonlinear Hamiltonian systems as well as separable ones, and even to non-Hamiltonian chaotic systems.

The essence of applying FPA to long-term computation is divided into two parts. One is the period detection of a unit cycle. The other is the computation of residue time equal to t − kT. This procedure helps to reduce the computation time for the long-term reliable simulation from O(t^1+1/M) to O(ln t + T / h₀). The FPA procedure improves not only the PMT method, but also facilitates the traditional symplectic method. The main problem of the symplectic method is that if the order M is not large enough (for example M<10) it still requires many computation loops for t = 10⁶⁰ –about O(60 ln 10 + T10^60/M) loops. For a medium-term time period, such as the ~10⁹ orbits of the Earth–solar system, the solution is required at a t = 10¹⁷ seconds magnitude. In this case, the symplectic method should work as well as FPA.

Recently, Barrio et. al[19] provided a shooting-periodical method which is a faster algorithm to obtain solutions in [0,10000] for a pre-obtained initial value for Lorenz system. The author notices that the long-term database obtained by Barrio is for the special initial values, while the method here is for any general initial values. It is well known that some Hamiltonian systems such as bounded Kepler system always have periodical orbits. Thus, in this case no shooting methods are needed to obtain the property initial values. While other Hamiltonian systems such as Arenstorf orbits[20] problems which have more freedom initial dimensions perhaps need such shooting methods to obtain precise initial values.

In this study, the author focuses on the long term simulation of periodic Hamiltonian systems, and not considers the spatial effect. When we count the spatial dimensional the ODE(ordinary differential equation) will change to PDE(partial differential equation), some of this PDE system still have periodic properties. For instant, the barotropic vorticity equation on the sphere[21]. While other systems may not have periodic for instant the Allee effect[22] in population dynamics. Two difficulties may occur when apply FPA to deal with these systems. Firstly, the determination of periodic of PDE system is more complicate then the ODE. Secondly, their may have no high enough integration scheme to do the integration within one periodic, and thus cause the founded periodic not very accurate. Nevertheless, as Sun et,al. pointed that the spatial dynamic pattern describe by PDE is ubiquitous in nature[23], thus improve FPA to do computation of such PDE dynamical system in an important work in the future study.

Acknowledgments

I thank Prof. Renguang Wu and Dr. Shuailei Yao for assisting with the evaluation of the paper.