Forward Period Analysis Method of the Periodic Hamiltonian System

Using the forward period analysis (FPA), we obtain the period of a Morse oscillator and mathematical pendulum system, with the accuracy of 100 significant digits. From these results, the long-term [0, 1060] (time unit) solutions, ranging from the Planck time to the age of the universe, are computed reliably and quickly with a parallel multiple-precision Taylor series (PMT) scheme. The application of FPA to periodic systems can greatly reduce the computation time of long-term reliable simulations. This scheme provides an efficient way to generate reference solutions, against which long-term simulations using other schemes can be tested.


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][2][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 9 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 −44 s) to the age of universe (about 10 17 s). Thus, a meaningful nondimensional time for a dynamical system is within 10 60 orders of magnitude. Simulation of the position of dynamical variables (not the trajectory structure) in the ultra long-term (i.e., t = 10 60 ) 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][8][9][10][11][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 The initial values are x 0 = 0 and p 0 ¼ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi , where x is the displacement and p is the momentum of the particles. Using the substitution y = e −x yields where y 0 = 1 and p 0 ¼ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1 À 0:02 p , and H ¼ p 2 2 þ 1 2 ðy 2 À 2yÞ is a constant. The Taylor series expansions relevant to solving (2) are where the coefficients are given by a k ¼ 1 and h is the step-size. The coefficients are determined from the initial conditions , and the relations 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]: , and g(p) = p. The second and more complex example is a non-separable Hamiltonian system, which is often broadly defined as: The Hamiltonian is H ¼ 1 2 p 2 þ pcos q, and the initial values are q 0 = 0, p 0 = 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: a mÀ i b i , d m = α m + g m , and the initial coefficients are α 0 = p 0 , β 0 = q 0 , b 0 = sin q 0 , g 0 = cos q 0 , c 0 = p 0 sin q 0 , and d 0 = p 0 + cos q 0 . The conservation of the Hamiltonian indicates dH dt ¼ 0. In the present numerical simulation, if |ΔH| 10 −16 (i.e., the smallest relative error in double-precision floating point arithmetic); it is regarded as unchanged. When H is a non-zero constant, j DH H j 10 À 16 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 where C 1 and C 2 are constants, and y N 2 2 À y 2 2 % jyðy N À yÞj jyjC 2 h Mþ1 . Since |p| and |y| are bounded variables, where C is a constant that satisfies |p|C 1 + |y|C 2 + C 2 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 highorder 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 7 , 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 5 time units. 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. 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.

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 −30 ) indicates the uncertainty (Δt) in estimating the period. Because E x is small and dxðtÞ , 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 −30 .
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, p ¼ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1 À 0:02 p and _ y ¼ À py ¼ À 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 2 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 −30 , h % 10 −15 must be set, and this requires 10 15 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 selfverification 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 Table 1. To the author's knowledge, this level of accuracy has never been reported.

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    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 30 − 2πk, to be more uncertain. Therefore, precision to 16 significant digits for sin(10 30 ) is not possible in a double-precision float platform.
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 −16 . The true residual time can be written as Under this situation, we have T r = 10 30 − 2πk, T 0 r ¼ 10 30 À 2pð1 þ εÞk and the difference between T r and T 0 r is Δt = 2πkε. The error bounds, ε, satisfying 2πkε < 10 −16 will guarantee Δx < 10 −16 . Since 0 < 10 30 − 2πk < 2π, 2πk % 10 30 and ε < 10 −16 /2πk = 10 −46 is the relative error bound of 2π. The 50 significant digits of π satisfy this error bound. Thus, T r = 3.231831977487846 and sin(10 30 ) = −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 30 ) 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 0 ) in the first stage, and to O(ln t) in the second stage, i.e., O(ln t + T / h 0 ).
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 0 ) at most. This improvement greatly decreases the CPU time cost and makes many unsolvable long-term problems be reliably solvable.
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 60 ) 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 5 are incorrect. It took about one day to obtain the result at t = 10 7 by direct integration, so obtaining a result at 10 60 is a seemingly impossible task for direct integration, but with the help of FPA, reliable results can be obtained. Another classical dynamical system is the mathematical pendulum. The Hamiltonian of a pendulum system is H ¼ 1 2 p 2 À cos q, and the dynamical equation is  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 −30 . This experiment again proves the correctness of FPA.
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 2 [0,10 60 ] 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 2 [0, 10 7 ], 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 0 ). 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 60 -about O(60 ln 10 + T10 60/M ) loops. For a medium-term time period, such as the~10 9 orbits of the Earth-solar system, the solution is required at a t = 10 17 seconds magnitude. In this case, the symplectic method should work as well as FPA.
Recently, Barrio et. al [19] provided a shooting-periodicalmethod 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.