Energy Minimization Scheme for Split Potential Systems Using Exponential Variational Integrators

: In previous works we developed a methodology of deriving variational integrators to provide numerical solutions of systems having oscillatory behavior. These schemes use exponential functions to approximate the intermediate conﬁgurations and velocities, which are then placed into the discrete Lagrangian function characterizing the physical system. We afterwards proved that, higher order schemes can be obtained through the corresponding discrete Euler–Lagrange equations and the deﬁnition of a weighted sum of “continuous intermediate Lagrangians” each of them evaluated at an intermediate time node. In the present article, we extend these methods so as to include Lagrangians of split potential systems, namely, to address cases when the potential function can be decomposed into several components. Rather than using many intermediate points for the complete Lagrangian, in this work we introduce different numbers of intermediate points, resulting within the context of various reliable quadrature rules, for the various potentials. Finally, we assess the accuracy, convergence and computational time of the proposed technique by testing and comparing them with well known standards.


Introduction
During the last few decades, there has been a growing interest in the development of numerical integration schemes for physical systems based on the theory of variations. Towards this aim, many authors suggested numerical methods that can be deduced through the discretization of the Hamilton's principle, leading to the class of integration methods known as discrete variational integrators. Due to their derivation process, they show specific geometric advantages, especially when they are applied to solving mechanical systems. It is worth noting that the preservation of the geometric features of the dynamical system makes these methods appropriate for both conservative and nearly dissipative systems [1][2][3][4].
When solving numerically ordinary differential equations (ODEs), one of the most difficult problems is the derivation of solutions appropriate for systems presenting highly oscillatory behavior [1]. For such systems, it has been shown that standard numerical methods require a huge number of (time) steps to track the oscillations, while geometric integrators proved to be quite advantageous tools [5]. The latter class of integrators endow highly qualified simulations for long term integration processes without showing spurious effects (i.e., bad energy behavior appeared in traditional schemes for conservative systems) [5][6][7][8][9][10].
On the other hand, in solving numerically Hamiltonian systems, an interesting still open issue is related with the property of the corresponding potential-energy function to be analyzed in several components each of them representing strongly varying dynamics. These interactions often (mainly) lead to extremely stiff potentials that force the solution to oscillate with a very small time scale [5]. To derive numerical integration schemes that respect the geometric characteristics of such systems, following our previous work [11], we split their potential energy into one fast and one slow component. For the resulting energy functional, we subsequently derive a special kind of high order geometric schemes based on exponential variational integrators. The latter schemes may employ different quadrature rules in the different potentials involved in the corresponding discrete action. For the sake of completeness of the proposed schemes, we assess their accuracy, convergence and computational time while comparing them with some standard ones.
The remainder of the article is organized as follows. We start by recalling the Hamiltonian formalism via the definition of the exponential integrators (Section 2). Then, we briefly summarize the physical concepts of the Lagrangian mechanics that help us to introduce the new exponential variational integrators (Section 3). These integrators are afterwards extended appropriately in order to be applicable to Hamiltonian systems of which the potential term of the energy functional can be split in fast and slow parts (Section 4). In Section 5, we present the tests and discuss the level of accuracy, convergence and computational time of the new schemes coming out of several numerical experiments performed on the general N-body problem. Finally (Section 6), we summarize the main conclusions extracted in the present work.

Exponential Integrators
In previous works we developed exponential integrators for Hamiltonian systems satisfying differential equations of the form [5] q + Ωq = g(q), where Ω denotes a diagonal matrix and U(q) is a smooth potential function. Due to the nature of the problems, the entries ω of Ω are assumed to have large modulus. In this work, we focus on the long time behavior of the numerical solutions, emphasizing on the special cases when ωh (with h being the time step) is relatively large [10]. Following Ref. [5], an exact discretisation of q(t) in (1) can be written as q n+1 − 2 cos(ωh)q n + q n−1 = 0, n ∈ Z Using the above anzatz, we can consider general numerical schemes defined as The latter expressions are fully determined via the functions ψ(ωh) and φ(ωh) which are assumed to be even and real-valued functions satisfying the condition ψ(0) = φ(0) = 1, see [5]. We note that, due to the presence of exponential terms in the above expressions, the resulting numerical schemes are known as exponential integrators (see, e.g., Appendix of Ref. [12]).

Variational Integrators Based on Interpolation Techniques
For the derivation of high order variational integrators we are interested in this work, we need to recall basic concepts from the discrete variational calculus [4,8,10]. Then, we define a discrete Lagrangian L d as a function of positions only, i.e., (Q represents a given configuration space). The discrete Lagrangian L d , is normally obtained from a continuous Lagrangian L(q,q) defined via the known map (T represents the tangent space of the configuration space Q). Towards this aim, an approximation must be considered of the form [4,10] (t k and t k+1 denote adjacent time moments of the time discretization grid, with t k+1 − t k = h).
In a similar way, the action sum S d is defined via the map Obviously, the action sum S d must correspond to a set of discrete Lagrangians L d , connected through the use of a discrete trajectory γ d = (q 0 , . . . , q N ) as [10] Subsequently, if we denote the dual space of TQ by T * Q, the latter space consists of covectors α such that In general, any covector α can be decomposed as A special case of this type of covectors are the dL d (q 0 , q 1 ) which, by applying the above definition, can be decomposed similarly as [4,10] where D i L d indicates the derivative with respect to the i-th argument of L d . Furthermore, according to the discrete variational principle, the solutions of the discrete system of equations determined by the Lagrangian L d must extremize the action sum S d of Equation (8) for given fixed points q 0 and q N . Extremizing S d over all the intermediate points of the discrete trajectory γ d , we obtain the system of difference equations The latter equations, are usually called discrete Euler-Lagrange equations [4,10]. Guided by the definition of the conjugate momenta in the case of the continuous mechanics, at this point we consider the discrete conjugate momenta at the time steps k and k + 1 as which are known as "position-momentum forms" of variational integrators. The latter are used to obtain, at the first step, the (q 1 , p 1 ) when the initial condition (q 0 , p 0 ) are provided [4].

High Order Variational Integrators
In order to derive (more accurate) high order methods, we approximate the action integral along the curve segment bounded by the two endpoints q k and q k+1 by utilizing a discrete Lagrangian that depends only on the given end points. Under these circumstances, we introduce expressions describing the intermediate configurations q j k and the corresponding velocitiesq j k , for j = 0, . . . , S − 1, S ∈ N, at time moments t j k ∈ [t k , t k+1 ], written as , with the special cases C 0 k = 0, C S−1 k = 1, for the time step h ∈ R. As a consequence, the above expressions for the intermediate time moments provide correspondingly the q j k andq j k given by In the particular case of oscillatory problems [10], we can subsequently define the intermediate configurations of the form It is worth noting that, through the chosen trigonometric expressions, the oscillatory behavior of the solution is inserted, see [10,11]. In addition, in order to satisfy the continuity condition, the equations must be fulfilled. Proceeding further, for any different choice of interpolation employed, we may define the discrete Lagrangian L d by a weighted sum of the form [10,12] where, for maximal algebraic order, we have proved that where m = 0, 1, . . . , S − 1 and k = 0, 1, . . . , N − 1 see [11]. For the application of the above interpolation technique combined with the trigonometric expressions of Equation (16), the crucial parameter u must be introduced. In our previous works, the determination of u came out of the phase lag analysis of Refs. [13][14][15][16]. In the latter works, we had showed that, for problems involving a domain frequency ω, the parameter u could be taken as u = ωh.
Moreover, for the solution of orbital problems of the general N-body system, where no unique frequency can be given, the parameter u may be defined by estimating the frequency of the motion of any moving individual body. The frequency estimation is necessary in order to find for each body (i) the frequency at an initial time t 0 and (ii) the frequency at later time moments t k for k = 1, . . . , N − 1. Following [10], we estimate the frequency of any moving mass using the definition of the curvature where q i (t), i = 1, ..., N, denotes the trajectory of the i-th particle [10].
In the present work, we would follow a different procedure and we will, first, reformulate the proposed integrators and, then, we combine them with the exponential integrators of Ref. [5] (for more examples see [17][18][19][20]). By working this way, we can get expressions for the parameter u which, afterwards, will be assessed on several individual cases of the general N-body problem.
The reformulation of the high order variational integrators of the present work offers the advantage of solving problems related with the Hamiltonian system (1). Towards this purpose, the discrete Euler-Lagrange Equation (12) deduced from the discrete Lagrangian (18) that relies on the intermediate configurations of (15) can be written as where By combining the latter variational integrator with that of (3), one may obtain the expression which, working in a similar way, defines exponential variational integrators based on the functions Ψ and Φ.

Determination of the Key-Role Parameter u
The high order exponential variational integrators of Equation (23), can be further elaborated in order to deduce an expression for Λ(u, ω, h, S) appropriate for determining the key-role parameter u introduced by the interpolation of Equation (16). In this way, exponential variational integrators that utilize the configurations q j k and velocitiesq or, similarly Λ(u, ω, h, S) + 2 cos(ωh) = 0.
Moreover, the feasibility of the crucial parameter u may be tested through Equation (25) and the consideration, as a concrete example, of the typical Lagrangian which corresponds to the simple harmonic oscillator with frequency ω By exploiting the weighted sum of Equation (18) as a representation of the discrete Lagrangian L d (q k , q k+1 ), in terms of the exponential variational integrators of (18) it takes the form The behavior of the exponential variational integrator obtained as described above, has been examined in this work for different choices of u and the results are illustrated in Figure 1. In this figure we considered the oscillation frequency to be ω = 1, the time  Even though the choice of u can be further examined, here we will consider it as coming out of Equation (25). This is especially advantageous for the extensions of the proposed techniques when they employed to treat Lagrangians of split potential functions, i.e., systems where the potential can be decomposed into several components involving different dynamical properties.

Family of Exponential High Order Variational Integrators
As an extension of the above method, in this section we derive a "family" of exponential high order variational integrators, following Ref. [11]. Towards this end, we split the potential energy of the Lagrangian into fast V f and slow V s terms, i.e., we write where T(q) denotes the kinetic energy of the system written as (M is a constant mass matrix) [21][22][23]. In this work, rather than using different quadrature rules to approximate the contribution of each component of the potential into the action (as done in Ref. [21]), we will follow a more general approach and we will use different number of intermediate points for the V f and V s [12,24,25]. By assuming that S 1 and S 2 represent the number of the points for the slow and fast potential, respectively, here we will restrict ourselves to choices for which S 1 < S 2 (the choice S 1 = S 2 creates the exponential integrators of Section 3).
Before embarking to specific simulation experiments, we should note that, after performing the discretization in (29) and adopting different quadrature rules for the two different potentials, the resulted discrete Lagrangian can be cast in the form where the intermediate positions q(t j k ) and velocitiesq(t j k ) may be taken from Equations (15) and (16), while the weights w 1 and w 2 could be deduced by exploiting Equation (19).

Numerical Examples
As representative examples, in this section we consider numerical solutions of the general N-body problem. As it is well known, the gravitational motion of an N-body system is described by a Lagrangian function of the general form [5,6] where G stands for the gravitational constant. In the following subsections, we are going to study two special cases of the latter Lagrangian function: (i) the corresponding to the satellite solar system (where N = 3), and (ii) the corresponding to the complete solar system (where N = 11).

Gravitational Motion of the Satellite Solar System
The satellite solar system (also called modified solar system) consists of three bodies (N = 3): the main (big) body and the two relatively small ones (planets) moving under the influence of their mutual gravitational attraction. For the description of this kind of planar motion, we consider masses m 1 = 1, m 2 = m 3 = 10 −2 , initial configurations q 1 = (0, 0), q 2 = (1, 0), q 3 = (4, 0) and initial velocitiesq 1 = (0, 0),q 2 = (0, 1),q 3 = (0, 0.5). Under these circumstances, it has been shown that the trajectories of the two planets (with respect to the main body, the Sun) are nearly circular with periods close to 2π and 14π, respectively (for more details see Ref. [5]).
In order to illustrate the advanced features of the proposed method, we considered, as a next step, an additional (fourth) body a satellite of mass m 4 = 10 −4 which moved rapidly around the mass m 2 . In this way, we included in the initial system a body moving with significantly different velocity (compared to that of the planets). The Lagrange function that described the final (four-body) system was then written as where In order to proceed further, we chose for the satellite m 4 the initial position and velocity to be: q 4 = (1.01, 0) and p 4 = (0, 0), respectively. Figure 2 shows the comparison of the global errors in the (four-body) system's total energy at t = 1 for time-step choices: h ∈ {10 −4 , 10 −2 , 5 × 10 −2 }, see [10]. Towards obtaining these results, we applied the proposed exponential variational integrators of Section 4 with the use of different numbers of intermediate points for the slow potential and the fast one of (34), namely S 1 (for the slow) and S 2 (for the fast). In fact, in Figure 2, the graphs correspond to the choices S 1 = 3, S 2 = 3 (blue dashed line), S 1 = 3, S 2 = 5 (red line) and S 1 = 5, S 2 = 5 (green line).
From Figure 2, it is obvious that the the graphs derived for S 1 = S 2 = 3 (blue line) were simply the standard high order variational integrator schemes of Ref. [10]. For all step sizes tested, the results obtained by increasing S 2 (while keeping S 2 > S 1 ), showed that the energy errors became progressively smaller. For the sake of completeness we have, in addition, tested the proposed methods with respect to their computational time consumed to obtain the same position errors. These simulated experiments showed us an interesting result which is depicted in Figure 3. As one can see, when comparing the method using S 1 = 3, S 2 = 5 (red-solid line) with the one using S 1 = 5, S 2 = 5 (green-dashed line), the latter scheme was more time consuming.

The Complete Solar System (N = 11)
As a final and more complicated example, we further considered the numerical solution of the complete solar system, see [5]. This is a special case of the N-body problem for N = 11, which consists of the nine planets of our solar system orbiting around the Sun and the Moon which is orbiting around the Earth (by adding a satellite in one of the planets, the mathematical problem becomes, obviously, more complicated). For the sake of accuracy and completeness, all initial conditions (positions and velocities) were taken from Nasa's JPL data base.
The general N-body system can be described through the Lagrangian for G = 2.95912208286 × 10 −4 [5]. In order to test the proposed schemes of Section 4, we split the potential of the latter Lagrangian in two terms as Initially, we chose S 1 = 3 for the slow potential and S 2 = 5 for the fast one, and study the time evolution of the various energy-components of the system. In Figure 4 we plot the various terms of the Lagrangian (32), namely the potential (black dashed line) and the kinetic energy (blue dashed line) components. Additionally, we plot the calculated total energy (red line), which was expected to remain constant throughout the integration process. In this simulated experiment, the time step was chosen to be equal to 2 days (a relatively long time step was required to keep on truck fast orbiting objects, the Moon in our example). In order to assess the accuracy and stability of the method, all the results presented were obtained for 10 6 days (a rather long time). The good and stable behavior of the method is evident from Figure 4.  Finally, the global errors for the total energy of the system at t = 1000 for time steps h ∈ {0.5, 1, 2.5, 5} days were compared. In Figure 5 we show the results obtained by applying the above schemes choosing S 1 = 3, S 2 = 5 (red line), S 1 = 5, S 2 = 5 (black line) as well as those obtained with S 1 = S 2 = 3 (blue line). From all cases tested, it is apparent that the smaller energy errors found corresponded to choices for which S 2 > S 1 , a conclusion which is in agreement with the findings of the previous simulated experiment. It is worth noting that the above numerical examples demonstrated the same errors in energy by using S 2 > S 1 compared to S 1 = S 2 = 3 for approximately 15% longer time step. This implies that, the proposed schemes had excellent energy preserving properties compared to standard schemes.
Before closing, it should be stressed that the proposed method is based on the high order variational integrator schemes of Ref. [10,11] that have been properly deduced to address nonlinear phenomena and specifically stiff problems in which some physical observable is varying with exceptionally high-frequency. The resulted discrete equations for such systems, may be viewed as coming out of a finite-dimensional Hamiltonian with special symmetry properties. Subsequently, these discrete equations are time-integrated with the aid of the new variational integrators. The latter advanced schemes possess enhanced numerical stability properties in the nonlinear regime (see also Ref. [16,25] where the above exponential variational methods with at least one intermediate point have been proven to be unconditionally stable).

Conclusions
In the present work, we developed a special type of numerical integration schemes which are based on the discretization of the Hamilton's principle. They are appropriate to be applied in systems that are described via an energy functional in which its potential component may be separated into various parts with strongly different dynamics. Interesting simulation problems of this kind appear in N-body gravitational systems, where their mutual interactions can be written as mentioned before and may lead to extremely stiff potential terms that force the solution to oscillate on a very small time scale.
The proposed high order numerical integrators rely on exponential functions entering the discrete action and on the application of different quadrature rules to the different component potentials. In this way, the derived integrators respect the geometric characteristics of the system while at the same time show great accuracy, convergence and short computational time compared to standard schemes.
Funding: This research received no external funding.