On stabilization of energy for Hamiltonian systems

Abstract

We study post-stabilization of numerical integrators for an autonomous Hamiltonian system in the form $H=H_1+H_2$, where $H_1$ and $H_2$ are two constants of motion. It is better to stabilize the two energies $H_1$ and $H_2$, respectively. As a particular case, the effectiveness of post-stabilization by the total energy is equivalent to one by the two independent energies when the splitting Hamiltonian becomes isotropic.

Introduction

Two problems should be focused on in the course of a numerical integration. The stability of differential equations of motion should be checked. For instance, Lyapunov's instability of Keplerian trajectory [1] would result in rapid increase of numerical errors. On the other hand, the preservation of a certain invariant manifold should also be considered because errors of the numerical algorithm may yield a great bias from the invariant.

For the first problem depending on the instability of equations in the sense of Lyapunov, it can be solved by virtue of a suitable coordinate transformation accompanied with regularization [2]. For example, the unstable equations of a two-dimensional Keplerian motion can be reduced to stable ones through the transformation of Levi-Civita and regularization. It is the same with equations of the three-dimensional Keplerian problem in terms of KS-transformation [2]. Rauch and Holman [3] noted that the second order leap-frog symplectic integrator of Wisdom and Holman [4] applied to the integrable Stark system with a larger eccentricity displays a chaotic behavior until the regularization is adopted. However, there is not a universal method to find adaptive transformations and regularization. As to the second problem, it relates to a practical problem how to satisfy relevant constraints in constrained systems where there are more coordinates than there are degrees of freedom. In principle, the number of numerically integrated first-order differential equations is demanded to be equivalent to the dimension of the phase space. This is called a reduction of the order of differential equations of motion by use of constraints [5], which is regarded as the method for satisfying the constraints automatically. However, the reduction is not valid in many cases because of the difficulty in choosing the sign after an operation of square roots from these constraints. As a particular case, a symplectic integrator can correspond to a modified integral of energy, but it cannot maintain other integrals of motion in general. In addition, it cannot be applied to non-Hamiltonian systems. Without loss of generality, it is a preferable and powerful technique for achieving the above two aims to employ Baumgarte's method [2], [6] called stabilization of differential equations. As one of stabilization methods for numerical integration, the post-stabilization [7], [8], [9] applies a stabilization step about the invariant at the end of each time step. It is usually superior to Baumgarte's method, since it not only maintains the accuracy of the invariant in magnitude of $\mathrm{O}({\tau }^{2(\sigma +1)})$, where τ and σ represent the time step and order of a numerical method, respectively, but also avoids encountering difficulty in choosing the best stabilizing parameters.

In general, there exist several constants of motion for a given dynamical system. In this case, it is vital to choose the best one of them as a stabilizing term because a stabilization method does not necessarily conserve all the constants. For an autonomous Hamiltonian system, stabilization by energy is interesting. The reason is that it is more important to guarantee the stability of some orbital elements including the semi-major axes, the mean motion, or the orbital period than to do the stability of other elements over a long-term integration interval. Some references [1], [5], [9], [10] noted the importance of energy correction, but we find also that Bartler [11] adopted the angular momentum as a constraint to estimate effects of ordinary differential equation stabilization on Lyapunov exponents. In view of these facts, a primary contribution of the present paper is to consider numerical post-stabilization of a Hamiltonian system with the separable form $H=H_1+H_2$, where $H_1$ and $H_2$ are all constants. In detail, we will trace whether the efficiency of post-stabilization by the total energy H or that by the two independent energies $H_1$ and $H_2$ is better.

The paper is organized as follows. Following Baumgarte [6] and Chin [7], in Section 2 we show how to operate the stabilization and post-stabilization. In this paper, the latter is particularly worth concerning owing to the convenience of application. To see its superiority, we use several simple examples in Section 3. First we apply the second order refined Euler scheme with and without post-stabilization to a two-dimensional linear oscillator with two constant frequencies. Here we are mainly interested in comparing the validity of post-stabilization between by the total energy and by the two independent energies when the system is isotropic (i.e. two identical frequencies) or anisotropic (namely two distinct frequencies). It will be found that there are some differences in the effectiveness of the two post-stabilization methods for anisotropic systems. To examine the fact further, then we apply other integrators such as a 4th order Runge–Kutta method to the linear oscillator. In particular, a detailed analysis of these differences is given. In addition, taking a planar Kepler problem and two-dimensional nonlinear coupled oscillator as representatives of isotropic and anisotropic systems, respectively, we test the effectiveness of the two aspects of post-stabilization again. Of course, it is necessary to employ an appropriate coordinate transformation such that a new Hamiltonian can be split into two constant pieces. Other similar physical models are listed also. Finally, we present our conclusion and discussion in Section 4.