1 Introduction

Constrained Hamiltonian systems have a lot of applications^[1-4],such as multi-body dynamics,classical mechanics,electrical circuits,molecular dynamics,chemistry,and fluid dynamics. The dynamic equations of constrained Hamiltonian systems are usually differential algebraic equations (DAEs). A lot of numerical methods can be used to solve DAEs,such as classical Runge-Kutta (R-K) methods,linear multistep methods,and generalized-α methods. However,the properties of the Hamiltonian system,such as symplecticity,energy,constraints,cannot be preserved. Therefore,more and more scholars have studied structure-preserving algorithms.

For DAEs,the drifts of the computed solution from the position,velocity or acceleration constraint manifolds are likely to occur during the simulation^[5-6]. One of popular techniques to overcome the difficulty is the constraint violation elimination technique whose main idea is the projection^{[5, 7]}. The projection approach has been proven to be valuable for a huge class of problems related to DAEs^[8]. Ascher and Petzold^[9] first proposed projected implicit RK methods for the solution of index-2 Hessenberg systems of initial and boundary value of DAEs and proved existence,stability,and superconvergence for the corresponding numerical approximations. Then,Ascher and Petzold^[10] extended these methods to higher-order higherindex Hessenberg systems and proposed a projected collocation method for the higher-order ordinary differential equations (ODEs) part of the DAEs and a projected invariant method for higher-index problems. Schropp^[11-12] proposed projected R-K methods for the solution of index-3 DAEs in the Hessenberg form and presented the corresponding convergence theorem,and he analyzed the geometric properties of these methods in the vicinity of equilibria,periodic orbits or asymptotically stable invariant sets. Chan et al.^[13] developed post-projected R-K methods,and the numerical approximation is projected only as part of the output process for index-2 DAEs. Compared with standard projections,these methods are strictly stable at infinity,and the order of convergence is unaffected.

The above mathematical studies for DAEs are mostly related to the constraint violation of the available integration schemes. Considering that the energy cannot be preserved by the previous projections,Hairer^[14] developed symmetric projection methods which were used to enforce conservation of energy for highly oscillatory Hamiltonian systems by Ascher and Reich^[15]. Andersen^[16] proposed the famous rattle algorithm which is symmetric,symplectic,and nearly preserves the Hamiltonian for molecular dynamics calculations by adding a projection step^[17]. But it is only of order two and thus not efficient for high accuracy requirements. Thus Console et al.^[17] presented symmetric linear multistep projection methods which extended the excellent rattle algorithm of order two to an algorithm of arbitrarily high order. They proved that symmetric linear multistep projection methods have the same qualitative behavior as the rattle algorithm. Calvo et al.^[18] used R-K methods and a projection technique to preserve any given Lyapunov function for ODEs. Calvo et al.^[19] put forward embedded projection methods by combining R-K formulae with a suitable projection. The aim of this projection was to preserve some first integrals in the numerical integration. Laburta et al.^[20] studied the numerical integration of non-conservative perturbations of differential systems by combining standard numerical integration methods with certain projection techniques. Hairer and Lubich^[21] put forward an energy-preserving R-K method for index-1 or index-2 DAEs by adding an energypreserving projection. Terze et al.^[22] recently proposed a Lie-group integration method for constrained multi-body systems in the state space. The mathematical model of multi-body systems dynamics is shaped as index-1 DAEs. In order to eliminate numerical constraint violations for generalized positions and velocities during the integration procedure,they introduced a constraint stabilization projection method based on the constrained least-square minimization algorithm.

The above projection methods are mostly used to study the problems of violations of constraints at the position and velocity levels or conservative of energy. However,the violation of constraint at the acceleration level has been neglected. Based on the augmented Lagrangian formulation (ALF),Bayo and Ledesma^[23] proposed mass-orthogonal projection methods for holonomic systems,which yielded no constraint errors at three levels (i.e.,at the position,velocity,and acceleration levels). Tian et al.^[24] extended the ALF with the mass-orthogonal projection method to the flexible multi-body dynamics based on the absolute nodal coordinate formulation. Dopico et al.^[25] extended the ALF to nonholonomic systems and discussed the effect of the velocity and acceleration projections on the constraint reactions and their evaluation. Malczyk and Fr¸azek^[26] presented a new parallel algorithm for the dynamics simulation of general multi-body systems. The ALF with mass-orthogonal projections was used to prevent from constraint violation errors. For the orthogonal projection,GarcÍa-Orden^[27] analyzed the relationship between the coordinate projection on velocities and the mechanical energy balance. The results of this analysis provided a practical criterion for the matrix selection. By combining the modified generalized-α method with the projection,Ding and Pan^[28] developed the generalized-α projection method which preserved constraints at three levels and the energy constraint. Although this numerical method maintains a high numerical accuracy for constraints at three levels,the energy error would have drifts.

Motivated by the advantages of aforementioned projected methods,in the present article,a projected R-K method,which preserves constraints at three levels and energy by combining R-K method with the projection technique,is developed. This paper is organized as follows. In Section 2,the equations of motion,which are DAEs,for constrained Hamiltonian systems are derived. In Section 3,a projected R-K method is established to solve DAEs. The validity of the proposed method by numerical results is presented in Section 4,and the concluding remarks are given in Section 5.

2 Governing equations of constrained Hamiltonian systems

The constrained Hamiltonian systems are modelled by n generalized coordinates q ∈ ℝ ⁿ subjected to l holonomic constraints,i.e.,

which are termed constraints at the position level^[29]. The Lagrangian function can be written as

where U (q) is the potential energy of the constrained Hamiltonian systems,λ ∈ ℝ^l is an array of Lagrangian multipliers,and T () is the kinetic energy and is considered to be a quadratic form of ,i.e.,

Note that in Eq. (3),M ∈ ℝ^n×n is an n × n symmetric positive-definite mass matrix. By using variational principles,the second Lagrangian equations can be written as

where q_i is the ith component of q.

Substituting Eqs. (2) and (3) into Eq. (4) yields the equations of motion for constrained Hamiltonian systems^[17],i.e.,

where ∇U (q) is the derivative of the potential energy U (q) with respect to q,and g_q(q) = is a Jacobian matrix of the constraint equations. By introducing a new variable of the velocity v = ,Eq. (5) can be transformed into first-order DAEs. In addition,for the sake of projection,the acceleration variable w = is introduced,and Eq. (5) can be rearranged as

For Eq. (6),there are many structure-preserving algorithms,such as symplectic R-K methods, Lie group methods,the generalized-α method,variational integrators,and combinations of these methods^{[2, 30-33]}. However,the constraints at three levels and the energy cannot be simultaneously preserved using these methods,which leads to invalidations of numerical simulations. To preserve constraints at three levels and energy,a projected R-K method will be presented in Section 3. Before that,some concepts must be introduced. By differentiating Eq. (1) once and twice with respect to the time,the constraints at the velocity level

and the constraints at the acceleration level

can be defined^[29]. Equations (7) and (8) should be conserved since Eq. (1) is held at any time. Another invariant for Hamiltonian systems is the total energy,

The integration results of Eq. (6) should also satisfy Eqs. (7),(8),and (9) theoretically. However,due to the unavoidable integration errors,the constraints at three levels and the conservation of total energy of the system would be violated during calculation. Therefore,projections have to be implemented to eliminate or minimize these violations.

3 Projected R-K methods 3.1 R-K methods

For unconstrained systems,R-K methods can be used directly to solve the first-order ODEs. Whereas,when the system has constraints,not all the state variables are independent^[34]. In this case,before utilizing R-K methods,the DAEs of the system must be transferred into firstorder ODEs. However,it is tedious to do this. For index-3 DAEs (such as Eq. (6)),the formulae of R-K methods can be found in Ref. [12],and Eq. (6) can be solved by the following s-stage R-K methods accordingly:

and Q_i,V_i,W_i ∈ ℝⁿ,Λ_i ∈ ℝ^l satisfy the nonlinear equations as follows:

where i = 1,2,· · · ,s. h is the time step,(∧) indicates that the solutions are obtained without the projection,b^T = (b₁,b₂,· · · ,b_s)^T ∈ ℝ^s and A = (aij)_s×s ∈ ℝ^s×s are the weight vector and the coefficient matrix of R-K method,respectively,and d^T = (d₁,d₂,· · · ,d_s)^T ∈ ℝ^s is defined as

If the projection is not performed,we just simply let

The coefficients A and b of R-K methods with the Butcher tableau can be written as

An s-stag R-K method of the Gauss type is symplectic and possesses the highest order 2s based on the Gaussian nodes for ODEs^[35]. Thus,a 3-stag R-K method of the Gauss type is adopted. The coefficients A,b,and c with the Butcher tableau can be written as^[33]

Note that when applying the algorithm (10) to Eq. (6),the nonlinear equation (11) must be solved at every step. Many methods are available for nonlinear equations,such as the Newton-Raphson iteration. Considering that the numerical solution of the algorithm (10) does not satisfy the constraints,which leads to a poor convergence^[36],we propose the following projected R-K methods in the next subsection.

3.2 Projected R-K methods

Numerical results obtained by the algorithm (10) do not preserve constraints at three levels and the energy. By using the projection technique,a new set of numerical solutions of position, velocity,and acceleration can be obtained,and these solutions satisfy constraints at three levels and conservation of total energy.

3.2.1 Orthogonal projection to constraint at position level

The position solution _n+1 does not completely satisfy the constraints at the position level,i.e.,Eq. (1). To satisfy Eq. (1) at t_n+1,we perform an orthogonal projection of the solution to the position constraint manifold and get q_n+1. It becomes the following constrained optimization problem:

The problem proposed by Eq. (16) can be solved by introducing the Lagrangian multiplier μ = (μ₁,μ₂,· · · ,μ_l)^T,and the Lagrangian function is

In order to minimize Eq. (16) according to the theory of constrained optimization,we differentiate Φ_n+1 with respect to q_ⁿ+1 and equate to zero,which yields

Combining Eq. (18) with Eq. (16),projection formulas can be achieved as

Solving the nonlinear equation (19),the projected position solution q_n+1 is obtained,and it satisfies the constraint at the position level exactly.

3.2.2 Orthogonal projection to constraint at velocity level

Similarly,the velocity solution _n+1 does not completely satisfy the constraint at the velocity level,i.e.,Eq. (7) at t_n+1. The orthogonal projection of solution _n+1 to the velocity constraint manifold is equivalent to the following constrained optimization problem:

where q_n+1 can be obtained by Eq. (13) or Eq. (19) depending on whether the projection of _n+1 is necessary. The Lagrangian multiplier ξ = (ξ₁,ξ₂,· · · ,ξ_l)^T is introduced to deal with this optimization problem,and the corresponding Lagrangian function is

Differentiating Ψ_n+1 with respect to v_n+1 and equating to zero yields

Combining Eq. (22) with the constraint equation (20),the projection formulas of velocity can be expressed as

3.2.3 Orthogonal projection to constraint at acceleration level

The same procedure can be applied to the acceleration solution _n+1. The orthogonal projection of _n+1 to the acceleration constraint manifold is equivalent to the following constrained optimization problem:

where G(q) is defined in Eq. (7). By introducing an array of Lagrangian multiplier ζ = (ζ₁,ζ₂,· · · ,ζ_l)^T,projection formulas of acceleration can be achieved as

3.2.4 Orthogonal projection to total energy

The projection to the total energy is a little different since there does not exit an energy state variable. The violation of energy can be attributed to the integration error of positions and velocities. Therefore,positions and velocities should be corrected so that the total energy can be preserved. By using the same procedure and introducing the Lagrangian multiplier η,the projection formulas to the total energy can be derived. However,there are two different energy projection strategies. The first one is widely used in the previous papers^{[21, 28]},which can be written as

The purpose of this correction is to equate H (q_n+1,v_n+1) with H (q_n,v_n). Whereas,the problem is that H (q_n,v_n) is not equal to H (q₀,v₀) due to the round-off error. As a result,the error of total energy will accumulate and increase with the time. Therefore,we propose another energy projection strategy,i.e.,

3.2.5 Multi-projection formulas

The formulas of the position projection,the velocity projection,the acceleration projection,and the energy projection have been given independently. However,all of these projections can be used at the same time.

If the position projection and velocity projection are used together,as most of previous literatures do^{[5, 9-13]},the formulas can be written as

And if all the projections are used,the formulas are

The only difference of Eqs. (29) and (30) is the last equation that whether equate H(q_n+1,v_n+1) with H (q_n,v_n) or with H (q₀,v₀). They are theoretically equivalent. However,they show different capabilities to preserve the total energy of the constrained Hamiltonian system numerically. Three numerical examples are used to explain the effectiveness of the proposed method.

4 Numerical examples

Example 1 Simple pendulum

In this section,a simple pendulum model composed of a point A with the mass m = 1 kg and a massless and nondeformable rod with the length l = 1 m is considered^[37]. And it is assumed that the mass point A is moving in the xy-plane as shown in Fig. 1. The generalized coordinates of the system are defined by q = (q₁,q₂)^T = (x,y)^T. Thus,the constraint equation can be written as g (q₁,q₂) = q₁² +q₂² −1 = 0. The gravity acceleration is assumed to be 10 m/s2 so that the potential energy is U(q) = 10_q2.

In order to make a comparison,four different methods are adopted to solve the equations of the motion of the simple pendulum system. The initial conditions are (q_1,0, q_2,0,v_1,0,v_2,0)^T = (0.8,−0.6,0,0)^T,and the time step is chosen as 0.001 s. Coefficients of R-K method can be seen in Eq. (15). The constraint Hamiltonian systems can be solved by projected R-K methods through the following procedures:

(i) Given the kinetic energy,the potential energy,and the constraints equations of the systems,the variational principle can be used to obtain the equations of motion of the systems. Furthermore,through the introduction of two new variables,i.e.,the velocity and the acceleration,respectively,Eq. (5) can be rewritten in the form of Eq. (6).

(ii) Solving the nonlinear equation (11) by the Newton-Raphson iteration with specific initial conditions,the time step and coefficients of R-K method in Eq. (15),we can get Q_i,V_i,W_i,Λ_i (i= 1,2,· · · ,s). Substituting these variables into Eq. (10),_n+1,_n+1,_n+1 can be calculated.

(iii) Four methods can be used to obtain q_n+1,v_n+1,and w_n+1,shown in Table 1,where R-K-PV represents the R-K method combined with the projections of position and velocity,R-K-PVAE-n represents the R-K method combined with the projections of position,velocity,acceleration and energy with respect to t_n,and R-K-PVAE-0 represents the R-K method combined with the projections of position,velocity,acceleration and energy with respect to t₀.

It is not so easy,or even impossible,to give analytical solutions for constrained Hamiltonian systems. Therefore,an important way to verify the numerical solution is to check whether the constraints and the total energy of the system can be preserved. The computational results of constraint errors at the position level,constraint errors at the velocity level,constraint errors at the acceleration level and total energy errors,are illustrated in Figs. 2-5,respectively,which are defined in the following equations:

Figure 2 presents constraint errors at the position level of four methods given in Table 1. It can be found that the constraint error at the position level of the R-K method dramatically increases after 6 s. It means the length of the rod varies from 0.6 m to 1.4 m with the time,which is apparently contradictory. The violation of constraint at the position level results in the invalidation of R-K method. In comparison,the errors of the other three methods maintain at the order of magnitude of 10⁻¹⁶.

The constraint errors at the velocity level by methods in Table 1 are plotted in Fig. 3. It shows that the error at the velocity level of R-K method increases steadily after about 6 s and arrives at the order of magnitude of 10⁸ at 30 s. It is much larger than other three methods which have both position projection and velocity projection.

The constraint errors at the acceleration level are depicted in Fig. 4. It is clear that the numerical solutions without the acceleration projection,i.e.,numerical solutions of R-K method and R-K-PV method,violate the constraints at the acceleration level. However,the violation of R-K method is much severer compared with the R-K-PV method which has position and velocity projections. The constraint errors at the acceleration level of R-K-PVAE-n method and R-K-PVAE-0 method remain the magnitude of 10⁻¹⁴.

The total energy errors of different methods are presented in Fig. 5. It can be seen that the total energy error of R-K method irrationally reaches the magnitude of 10¹⁵ after 30 s and keeps increasing super linearly,and the error of R-K-PV method is of 10⁻⁶ and increases linearly. In contrast,the R-K-PVAE-n method and R-K-PVAE-0 method,by taking the energy projection into account,show excellent performance in preserving total energy of constrained Hamiltonian systems. However,the R-K-PVAE-0 method is more accurate and steady in terms of preserving total energy in a long period of time.

Combining Figs. 2-5,it can be found that the R-K method without the projection fails to predict the motion of constrained Hamiltonian systems in a short time due to the violations of constraints at every level and non-preserved energy. By considering the position projection and the velocity projection,the R-K method gives a better result of constrained Hamiltonian systems,but still produces violations of constraints at the acceleration level and the dissipation of total energy. The R-K method with the position projection,the velocity projection,the acceleration projection and the total energy projection performs the best capability to preserve constraints at every level and total energy of constrained Hamiltonian systems. However,the projection method in Eq. (30) is more preferable.

In order to show the broad applicability of the present approach,two more numerical examples have been studied in the following manuscript.

Example 2 Double spherical pendulum

In Fig. 6,the double spherical pendulum is shown^[38-39],where l₁ and l₂ are the lengths of massless rods,and the masses of points A and B are m₁ and m₂,respectively. The generalized coordinates of the double spherical pendulum are introduced by q = (q₁^T,q₂^T)^T,where q₁ = (x₁,y₁,z₁)^T and q₂ = (x₂,y₂,z₂)^T are coordinates of mass points A and B. Similar to the example 1,the nonlinear constraint equations are defined as follows:

The total energy of the double spherical pendulum is calculated by

By using the proposed method,relative constraint errors of the double spherical pendulum at the position level,relative constraint errors at the velocity level,relative constraint errors at the acceleration level,and the relative energy error are,respectively,defined in the following equations:

where i = 1,2.

In this numerical example,l₁ = l₂ = 1 m,m₁ = 1 kg,m₂ = 2 kg,the gravity acceleration is assumed to be 10 m/s²,the vectors of initial coordinates q₀ = and initial velocities v₀ = (0,0,0,1,1,2) are taken.

Figure 7 shows the relative trajectory of the mass point B to the mass point A. From these figures,compared with the results of R-K,the relative trajectory obtained by the proposed method is well preserved. This shows the present method has the advantage of structurepreserving property.

Meanwhile,the relative constraint errors at the position level,the velocity level,and the acceleration level of the massless rods of l₂ are plotted in Figs. 8,9,and 10,respectively. It is found that the amplitudes of these errors obtained by the proposed method are very small. One needs to highlight that the proposed method not only can make the constraint errors at the position level and the velocity level but also at the acceleration level more well preserved,which can be seen in Fig. 10. What is more,the relative energy error of the double spherical pendulum is shown in Fig. 11. One can clearly find that the amplitude of the result obtained by the R-K-PVAE-0 is smaller than the other three kinds of methods which are mentioned in this manuscript. This shows that the proposed method can maintain not only the constraint errors but also the energy error well preserved.

Through this example,one can find that the proposed method can be also applied to more multiple degrees of freedom of complex system.

Example 3 Closed-loop rotating triangular tethered satellite system

In order to highlight the effectiveness of the proposed method and the broad applicability of the present approach,another very interesting study is carried out to the closed-loop rotating triangular tethered satellite system.

In Fig. 12,a closed-loop rotating triangular tethered satellites system is shown^[40]. There are some assumptions as follows: (i) the satellites are regarded as mass points; (ii) tethers are considered as inextensible,tight,and massless; (iii) the closed-loop rotating triangular tethered satellite system is rotating.

The vector of the generalized coordinates is defined by q = (q₁^T,q₂^T,q₃^T)^T,where q₁ = (x₁,y₁,z₁)^T,q₂ = (x₂,y₂,z₂)^T,and q₃ = (x₃,y₃,z₃)^T. Similar to the examples 1 and 2,the constraint equations can be written as

The total energy of the system is calculated by

where μ = 3.986 × 10¹⁴ m is the gravitational parameter of the Earth.

In this numerical example,the lengths of the tethers l_i (i = 1,2,3) are 1 km,and the weights of satellites m_i (i = 1,2,3) are 100 kg. Meanwhile,the initial coordinates and initial velocities of the system’s center of mass (CM) are (r₀,0,0) and ,where r₀ = 6.7 × 10⁶ m is the orbital amplitude.

The planes of m₁,m₂,and m₃ are parallel to the yz-plane in the initial time. The system is assumed that the initial rotation rate about the CM is .

This example is different from the examples 1 and 2 that the nonlinearity of the gravity is considered,which can be found in the last terms of Eq. (42). The three satellites have similar situations of the three constraints. For the sake of the brevity,the constraint of m₁ and m₂ is used as an example in the following manuscript.

The trajectory of the CM is plotted in Fig. 13. One can find that the trajectory is still stable,and no drifting. This can be proved that the proposed method has the advantage of maintaining the dynamic system stability.

Figures 14-16 show the relative constraint errors at three levels. From these figures,it can be found that relative constraint errors of l₂ at three levels remain in the small magnitudes of 10⁻¹⁵,10⁻¹¹,and 10⁻¹³,respectively. Furthermore,one can notice that as the time dramatically increases,the amplitudes of these errors still keep in a small range,which proves that the present method has the advantage of the long term stability. In Fig. 17,the relative energy error of the closed-loop rotating triangular tethered satellite system has been plotted. From Fig. 17,it can be observed that the amplitude is also small,which leads to the conclusion that the effectiveness of present method is not only suitable for the linear case (the examples 1 and 2),but also for the nonlinear case (the example 3).

5 Conclusions

By combining the R-K method with the projection technique,a method that preserves constraints at three levels (at the position level,the velocity level,and the acceleration level) and the total energy for constrained Hamiltonian systems is proposed in this paper. A 3-stage,6th-order implicit R-K method is adopted to solve the DAEs derived by the standard Lagrangian multipliers method,and the projection is then applied to correct the answer to satisfy the constraint equations at three levels and to preserve the energy. By numerical analysis,the proposed method exhibits good capabilities to preserve constraints not only at the position and velocity level,but also at the acceleration level and the total energy of constrained Hamiltonian systems. Meanwhile,the numerical experiments show that the advantage of the proposed method is long term stability. Furthermore,it is also found that the presented method can be extended to the nonlinear case and more complex problems.