1 Introduction

Many numerical methods have been developed to simulate the dynamics. Typically,the symplectic methods^[1-2] preserve the symplectic structure of the original flow of the differential system and have an excellent long-time energy behavior. During the last decades,they have been successfully applied to problems in mechanics,physics,biology,chemistry,etc.

There are two ways to construct symplectic methods. One is that the construction of symplectic methods is based on dual differential equations,the resulting methods,e.g. the precise integration methods^[3-4],symplectic Runge-Kutta (SRK) methods^[5-6],symplectic partitioned Runge-Kutta (SPRK) methods^[7-8],symplectic Runge-Kutta NystrÖm (SRKN) methods^[9-10],symplectic exponentially fitted Runge-Kutta (SEFRK) methods^[11-12],which have been developed for a large class of problems. Another way is symplectic methods derived from the discrete variational formulation of a given system. The resulting numerical integrators,known as variational integrators^[13],have remarkable characters,such as symplectic and momentumpreserving,as well as the superior predictions about energy conservation or change. Based on different discrete variational formulations (e.g. Hamilton’s variational principle,Lagrange d’Alembert principle,Hamilton-Pontryagin principle),variational integrators have been developed for conservative mechanical systems^[14-15],constrained systems^[16-20],and stochastic systems^[21-22]. The applicability of variational integrators is not limited to mechanical systems,but extended to the simulation of heat conduction^[23],electric circuits^[24],and magnetic systems^[25].

The theories of discrete Lagrangian mechanics and the corresponding variational integrators are reasonably well developed. But for some systems,for example,degenerate Lagrangian systems,the variational integrators may be more easily derived from the Hamiltonian side. However,there are not sufficient corresponding researches^[26-27]. And the Hamilton’s variational principle in the phase space corresponding to the boundary value problems is not convenient to construct variational integrators directly. Luo et al.^[28] proposed an unconventional Hamilton’s variational principle that absolutely corresponds to initial value problems and the resulting algorithms for the liner dynamic system.

In this paper,we propose a new approach for constructing variational integrators for conservative mechanical systems based on a simplified unconventional Hamilton’s variational principle in the phase space. And some particular variational integrators are derived. Meanwhile,the symplecticity and the computational performance of the proposed integrators compared with those of the SRK methods are discussed. Finally,we carry out some numerical examples to show the numerical behavior and the efficiency of the proposed integrators when they are compared with those of the SRK and the fourth-order Runge-Kutta (RK4) methods.

2 Simplified unconventional Hamilton’s variational principle

In the time step [0,τ],the functional of the unconventional Hamilton’s variational principle for Hamiltonian systems is as follows:

where q and p are the d-dimensional generalized displacement and momentum,respectively,H is the Hamilton function,q₀ and p₀ are the initial value vectors,q₀ = q(0),p₀ = p(0),q_τ = q(τ),p_τ = p(τ),the variables with (◦) are constrained variational variables which are regarded as invariants during the variational operation,and (·) indicates the differentiation with respect to the time t.

Using δⅡ(q,p) = 0,the Hamilton dual equations

and the initial conditions

can be derived,which clearly indicate that the unconventional Hamilton’s variational principle can reflect all the features of the initial value problems.

If q and p satisfy the initial conditions,Eq.(1) can be simplified as

From ,we obtain

Combined with the integrand by parts,we have

Considering the arbitrariness of δq and δp,the Hamilton dual equations can be easily derived. Equation (5) is the simplified unconventional Hamilton’s variational principle in the phase space presented in this paper,and Eq.(4) is the corresponding functional. q and p have satisfied the initial conditions. Therefore,Eq.(5) absolutely corresponds to the initial value problems.

3 Variational integrators based on simplified unconventional Hamilton’s variational principle 3.1 Construction of variational integrators

In the time step [0,τ],t₀,t₁,· · · ,t_m are a set of different specified points,and t₀ = 0. q and p are approximated by the mth-order Lagrange polynomial,which can be expressed as

where ,and q_i and p_i are the generalized displacement and momentum at the corresponding points,respectively.

Substituting Eq.(7) into Eq.(6) and using the Gaussian quadrature,the following relationships are established:

where g₁,g₂,· · · ,g_n are the Gaussian quadrature points,and b_j are the weight coefficient of Gaussian quadrature.

Considering δq₀ = 0,δp₀ = 0,and the arbitrariness of δq_k,δp_k,Eq.(8) can be rewritten as

q₀ and p₀ have already been known,therefore,there are 2md equations and 2md unknowns in Eq. (9). For linear or nonlinear dynamic systems,Eq. (9) is a set of linear or nonlinear algebraic equations.

Solving Eq. (9),the values of q_i,p_i (i = 1,2,· · · ,m) can be achieved. And the values of q_τ and p_τ can be derived from Eq. (7). However,the preferable method is introduced to solve q_τ and p_τ as follows. Computing the integral for Eq. (2) yields

Using the n point Gaussian quadrature results in

The values of q_τ and p_τ can be easily obtained from Eq. (11). It should be pointed out that Eq. (11) is also used to calculate the end point values in the SRK methods.

In this paper,the Newton-Raphson method is used to solve the nonlinear algebraic equations. Therefore,the selection of initial values is critical during the iteration process. Because of the particularity in constructing the integrators,the initial values for the first time step are arbitrary and for other time steps can be easily achieved by the extrapolation of the Lagrange interpolation polynomial of the last time step,namely,

where q_i and p_i are the values of q_i and p_i of the last time step.

The processing method can make the initial values much closer to the solution and decrease the iterations at the cost of slightly increasing the calculation of selecting initial values. And the numerical examples in Section 5 show that the processing method can make a 40%−50% efficiency improving.

3.2 Optimal variational integrators

In the following,the optimizing of the proposed integrators is discussed. Assuming m = n,the Lagrange interpolation points t₁,t₂,· · · ,t_m are consistent with the Gaussian quadrature points g₁,g₂,· · · ,g_n. Considering l_i(t_j) = δ_ij,δ is the Kronecker delta,and Eq. (9) can be simplified as

where

Solving Eq. (13),the values of q_i,p_i (i = 1,2,· · · ,m) can be obtained. And q_τ and p_τ are given by

which completes the calculation of a time step.

Equations (13) and (14) are the optimal forms of the proposed integrators in this paper. It can be seen clearly that the optimal form is much more concise and easier to program and for application.

4 Performance of proposed integrators 4.1 Symplecticity of proposed integrators

The method to prove the symplecticity of the proposed optimal integrators is described in the following. Letting v = [q^T,p^T]^T,F (v) = [H'_p^T,−H'_q^T]^T,Eqs. (13) and (14) can be,respectively,expressed as

Taking partial derivatives of Eqs. (15) and (16) with respect to v₀,we get

where given by Eq. (17) is substituted into Eq. (18),and is achieved. If A satisfies

where A is a symplectic matrix,and the proposed integrators are symplectic,,0_d and I_d are the d-dimensional zero matrix and unit matrix,respectively,and J_d is a standard symplectic matrix.

The above approach is the generalized way to prove the symplecticity of the proposed integrators. For different problems,the specific form of A can be obtained. It is shown that A is a strictly symplectic matrix from enormous examples. Therefore,the proposed integrators are symplectic.

4.2 Comparison of performance with SRK methods

In the following,we analyze the performance of proposed integrators and the SRK methods.

The m-stage 2m-order SRK method can be expressed as

The Jacobi matrix of the Newton-Raphson iteration method for the proposed integrators can be obtained from Eq. (15),which can be expressed as the combination of the linear part and the nonlinear part. The nonlinear part is as follows:

which is a block diagonal matrix,and the size of submatrix is 2d × 2d. However,for the SRK method,the nonlinear part of Jacobi matrix of Eq. (20) is a full matrix.

It can be seen that the proposed integrators compared with the same order SRK method have two advantages. One is that the nonlinear part of the Jacobi matrix of Newton-Raphson iteration method is a block diagonal matrix and that of the SRK method is a full matrix,which means less data that need to be updated during the iteration process and the less computational cost of the proposed integrators. Another advantage is that the block diagonal nonlinear part is much more helpful to construct more efficient iteration methods.

5 Numerical examples

In this section,several numerical examples are tested to access the numerical performance of the proposed integrators.

5.1 Linear harmonic oscillator

To study the symplecticity of the proposed integrators,the linear harmonic oscillator is considered. The Hamilton function is H = p²/2 + q²/2. For τ = 1,m = 2,the proposed integrators with 2,3,4,and infinite integration points are used,respectively,and the results are listed in Table 1.

The results show that,for initial value problems,the integrators based on the simplified unconventional Hamilton’s variational principle in the phase space are symplectic only under certain conditions,such as m = n in this paper. According to Table 1,the integrator with the infinite integration point in Ref. [28] is not strictly symplectic.

5.2 Perturbation pendulum problem

The perturbation pendulum problem with the Hamilton function H = p²/2−cos(q(1−p/6)) is an inseparable Hamiltonian system because of the coupling term. The initial conditions are q(0) = 1 and p(0) = 0.1. The simulations are carried out for 0 s ≤ t ≤ 10² s with time steps τ = 0.1 and τ = 0.2.

Firstly,we study the accuracy of the proposed variational integrators. According to Ref. [29],the accuracy order r of numerical algorithms can be evaluated as follows:

where τ₁ and τ₂ are the time steps,e₁ and e₂ are the corresponding errors of the integral results.

Table 2 shows the changes of accuracy order r with m,which indicates that the proposed integrators have accuracy of 2m order.

Then,the comparison of the accuracy between the proposed variational integrators and the SRK methods is discussed. In the cases of m = 1,2,3,and 4,the proposed integrators and the SRK methods are denoted as M2,M4,M6,M8 and SRK2,SRK4,SRK6,SRK8,respectively. Figure 1 illustrates the relative error e_H = |(H(t) − H(0))/H(0)| in the numerical solution of the Hamilton function versus the time t in the case of τ = 0.2. The figures show that the e_H curves of the two methods almost overlap for the same m and further indicate that the accuracy order of the proposed integrators is 2m,which is the same as the SRK methods.

Finally,the iteration efficiency of the proposed integrators with the initial values estimated by Eq. (12) is studied. The average iterations per step of the M4 and M6 integrators are shown in Table 3. The results indicate that the average iterations are much less,and the efficiency is improved by 40%−50%.

5.3 Classical trajectory of diatomic system

The classical trajectory method is an effective theoretical method in studying the microcosmic reaction dynamics. The micro-reaction system,as a classical Hamiltonian system,can be regarded as a particle system of which the atoms are moving at the electronic potential energy surface. Considering the diatomic system,the dimensionless Hamilton function is H = p²/2 + (e^−2q − 2e^−q)/2. The initial conditions are q(0) = 0,. The M4,SRK4,and RK4 methods are used for the calculation. The simulations are calculated to 10⁶ s with a time step τ = 0.1 s.

Figure 2 illustrates the relative error e_H with the time t. Obviously,the M4 and SRK4 methods with almost the same e_H curves both show the advantage of the long-time energy conservation,which accords with the property of symplectic algorithms. The e_H curve of the RK4 method increases linearly with t − t₀. Figure 3 describes the comparison of phase diagrams between the three methods. It is clearly shown that the calculating phase diagrams of the M4 and SRK4 methods are overlapped closed curves,which means that the two atoms vibrate periodically,and the numerical results are consistent with experiment results. The phase diagram of RK4 method is a closed loop with a certain width,which indicates that an obvious artificial dissipation phenomenon has appeared during simulation and also reflects the disadvantages of non-symplectic algorithms.

5.4 Kepler (two-body) problem

The celestial mechanics often involves the Kepler problem,for example,satellites circle the earth,and the planets move around the sun or binary star systems. The Hamilton function of the Kepler problem is The initial conditions are q₁(0) = 0.5,q₂(0) = 0, p₁(0) = 0,and ,and the period of solution is 2π. The M4,SRK4,and RK4 methods with the time step τ = 2π/256 are used to obtain the solution until t = 200π s. Letting z = [q₁,q₂,p₁,p₂]^T,the global error is measured by

The comparisons of the relative error e_H during the last 10 periods are shown in Fig. 4. The diagrams of global errors e_z versus the period T for various methods are illustrated in Fig. 5. It is obvious that the M4 and SRK4 methods both have the advantage of energy conservation in the long-time simulation,and the relative error e_z is stable of the order 10⁻⁷. The global error e_z of the M4 and SRK4 methods grows linear with t − t₀ and that of the RK4 method quickly grows with (t − t₀)² because of the energy dissipation. Therefore,symplectic algorithms,such as the M4 and SRK4 methods,show a huge advantage in the quantitative problems,especially in the long-term tracking numerical simulations.

5.5 Two-fixed-center problem

Consider the two-fixed-center problem. This problem with three degrees of freedom can be regarded as a simplified model for the restricted three-body problem model,which is of great importance in studying the movement theory of artificial earth satellite and not solved completely so far. The Hamilton function is ,where and μ satisfying 0 ≤ μ ≤ 1 is a parameter.

Take μ = 0.5,τ = 0.1. The initial conditions are q₁ = 1.0,q₂ = 0.1,q₃ = 0,p₁ = 0,p₂ = 0.5,and p₃ = 0.2. The problem is solved using the M4,SRK4,and RK4 methods,respectively,and carried out to 10⁴ s.

Figure 6 shows the relative error e_H with the time t. The figure indicates that,the M4 and SRK4 methods with a stable e_H of the order 10⁻⁶ have a preferable long-term energy conservation behavior,and e_H of the RK4 method is accumulated linearly along the time t− t₀.

6 Conclusions

In the present work,a new way for constructing variational integrators has been developed. The new integrators are based on the simplified unconventional Hamilton’s variational principle in the phase space that not only corresponds to initial value problems,but also is very convenient for application. For the nonlinear system,the selection of initial values in iterations for the Newton-Raphson method is proposed.

Numerical studies show that the variational integrators have 2m order computational accuracy,where m is the order of the Lagrange polynomial. Compared with the same order SRK methods,the computational cost of the proposed integrators is much less. In addition,the nonlinear part of Jacobi matrix has a block diagonal form,which helps to construct more efficient iterative algorithms.

Experimental tests show that the proposed integrators have the high accuracy and excellent long-time energy behavior. Meanwhile,it becomes obvious that the variational integrators based on the unconventional Hamilton’s variational principle are symplectic only if the interpolation points and numerical integration points have particular relations. Therefore,it is not entirely accurate that the algorithms constructed based on the unconventional Hamilton’s variational principle are all symplectic (see Ref. [28]).