- Email: [email protected]
A numerical method that conserves the Runge-Lenz vector
CHINESE ASTRONOMY AND ASTROPHYSICS
ELSEVIER Chinese Astronomy and Astrophysics 30 (2006) 87 101
A Numerical Method Runge-Lenz
that Conserves Vector t *
the
LIU Fu-yao 1,4 W U Xin 2'3 LU Ben-kui 1'4 1purple Mountain Observatory, Chinese Academy of Sciences, Nanjing 210008 2School of Science, Nanchang University, Nanchang 330047 3Department of Astronomy, Nanjing University, Nanjing 210093 4National Astronomical Observatories, Chinese Academy of Sciences, Beijing 100012
A b s t r a c t An exhaustive discussion is carried out on isolating integrals and the trapezoidal formula which can conserve the Runge-Lenz vector. An isolating integral is an invariant that restricts the region of particle motion. The autonomous integrable Hamiltonian system with n degrees of freedom has only n mutually involutive independent isolating integrals, and the existence of other isolating integrals is meaningful to the particle motion. In the Kepler two-body system there exist the energy integral, the angular momentum integral and the RungeLenz vector. These correspond to 3 independent isolating integrals in the case of plane motion, and to 5 in the case of space motion. In the former, the integrals makes up the symmetry group SO (3) of the system, which can be transformed into the symmetry group of the two-dimensional isotropic harmonic oscillator through the Levi-Civita transformation, which is accurately conserved by the trapezoidal formula. On the other hand, in the case of space motion, the strict conservation of the energy and angular momentum inegrals and the Runge-Lenz vector by the trapezoidal formula is manifested in the 5 Kepler orbital elements a, e, i, ft and w. K e y w o r d s : celestial mechanics--methods: numerical
t Supported by National Natural Science Foundation Received 2004 11 01; revised version 2004-11-22 * A translation of Acta Astron. Sin. Vol. 46, No. 3, pp. 294-306, 2005 0275-1062/06/S-see front matter © 2006 Elsevier B. V.All rights reserved. DOI: 10.1016/j.chinastron.2006.01.007
88
LIU Fu-yao et al. / Chinese Astronomy and Astrophysics 30 (2006) 87-101
1. I N T R O D U C T I O N In terms of linear and nonlinear relations in kinetic differential equations, Hamiltonian systems can be divided into two kinds, linear and nonlinear systems. The linear systems and some of the nonlinear systems belong to integrable systems, possessing the dynamic properties of ordered motion free of chaotic behaviour. The necessary and sufficient condition for an autonomous Hamiltonian system with n degrees of freedom to be integrable is that there exist n mutually involutive, independent isolating integrals [1,2]. Of course, some integrable Hamiltonian systems of n degrees of freedom possess more than n independent isolating integrals. For a linear system of n degrees of freedom with frequency commensuration, there exist 2n - 1 independent isolating integrals;-- one such example is the isotropic autonomous harmonic oscillator of n degrees of freedom. For the Kepler two-body system which pertains to an isotropic system, there exist the energy integral, the angular momentum integral and the Runge-Lenz vector. The first two respectively reflect invariance with respect to time translation and space rotation. In the case of plane motion, and all three represent independent isolating integrals, while in the case of space motion there are only 5 independent integrals [a] out of a total of 7 integrals (i.e., 1 of energy, 3 of the components of angular momentum and 3 of the components of the Runge-Lenz vector). A good numerical method should have a good computational efficiency, a high numerical accuracy as well as a good numerical stability and, especially, can conserve the geometrical topological structure and proper invariant manifold of the original system. Examples are the symplectic algorithm[ 4-7], the stabilization method and the method of reduction of order of equation Is,9]. Unquestionably, the often used symplectic algorithms such as the first order symplectic Euler method and the second order Wisdom-Holman leap-frog symplectic mapping, etc., are able to maintain the energy integral[ 1°], but this does not mean that the symplectic algorithm strictly conserves the energy integral. As a m a t t e r of fact, there always exists the truncation error In] when the Hamiltonian system obtained from the numerical results is compared with the original Hamiltonian system. However, the error in the Hamiltonian function has no long-term variation, so that the integration always proceeds near the plane of energy. Thus, the symplectic algorithm conserves an altered form of the energy integral, that is to say, the symplectic algorithm is a method which automatically satisfies the energy constraint condition. From this angle it is seen that the symplectic algorithm generally needs no stabilization treatment [12]. Besides, for the isotropic harmonic oscillator and the Kepler two-body system, Yoshida found that the symplectic algorithm can conserve the angular momentum integral in an altered form [13,14]. However, the symplectic algorithm, generally speaking, can not conserve all the isolating integrals. Yoshida pointed out that when the symplectic algorithm is applied to a linear anisotropic system with frequency commensuration, no other integrals in their altered forms exist except the energy integral in its altered form [131. He made further discussions and pointed out that the symplectic algorithm is unable to conserve the proper Runge-Lenz vector of the Kepler two-body system [14]. However, there is another sort of symplectic algorithms which can strictly conserve the energy integral and other integrals of motion. The simplest trapezoidal formula in the Newton-Coates family of integration as a linear multi-step implicit symplectic algorithm can accurately conserve the energy of the linear oscillatorI15-171. Here, it should be emphasized that the linear multi-step symplectic algorithm and the symmetrized algo-
LIU Fu-yao et al. / Chinese Astronomy and Astrophysics 30 (2006) 87-101
89
rithm are actually approximate symplectic structure preserving algorithms, i.e. subject to the conditions that the algorithm is stable and that the characteristic equation is determined by the principal root while the effect of spurious roots can be neglected [15-17]. For the linear Hamiltonian system, a sort of multilayer symplectic integrator in the shape of the trapezoidal formula can be constructed by making use of the Newton-Coates formula and the implicit midpoint formula [ls-2°]. As for the nonlinear system, Minesaki and Nakamura studied the feasibility of the application of the trapezoidal formula. Using the Levi-Civita or KS transformation, they transformed the nonlinear Kepler two-body system into the isotropic harmonic oscillator, to which they applied the trapezoidal formula. They then worked back to the original system, thus with the energy integral, the angular momentum integral and the Runge-Lenz vector are strictly conserved [21,221. Considering that the deduction by Minesaki and Nakamura of strict conservation of the Runge-Lenz vector with the trapezoidal formula is somewhat tedious, we seek in this paper a simpler derivation following a detailed exploration. The integrals of the Kepler two-body plane motion make up the symmetry group SO (3) of the system, which can be transformed into the symmetry group of the two-dimensional isotropic harmonic oscillator via the Levi-Civita transformation and the latter symmetric group can be strictly conserved by the trapezoidal formula. This way of derivation is comparatively simple. At the same time we shall make an in-depth exploration of isolating integral and the condition of integrability. We shall also examine the numerical accuracy of the six Kepler orbital elements. Five of the elements, i.e. a, e, i, ~ and ~, should be constants and we shall investigate whether they are strictly conserved by the trapezoidal formula. We shall also explore the condition under which the isolating integrals are conserved by numerical integration of the nonlinear integrable system including the two-dimensional coupled harmonic oscillator.
2. S Y M M E T R Y G R O U P OF T H E K E P L E R T W O - B O D Y M O T I O N For later convenience, we state in this section the concepts of isolating and non-isolating integrals as well as the symmetry group of the Kepler two-body plane motion, which can be transformed into that of the two-dimensional isotropic harmonic oscillator.
2.1 Isolating and non-isolating integrals Research on dynamical systems generally revolves around such concepts as isolating integral, non-isolating integral, integrability and non-integrability, etc. An isolating integral is the very invariant that restricts the region of motion of the particle. For an autonomous Hamiltonian system H = H ( q , p), if a function (I) = (I)(q, p) (of the generalized coordinate and momentum) is such that its Poisson bracket with H satisfies the involutory relation,
{e, H} = 0,
(1)
then (P is said to be an isolating integral of H. An autonomous Hamiltonian system of 3 degrees of freedom has at least one isolating integral because H itself is an invariant of motion (called the energy integral). It has at most 5 independent isolating integrals, because the number of independent variables in the phase space is five, and each isolating integral is a five-dimensional hypersurface. When the system possesses 3 or more independent isolating
90
LIU Fu-yao et al. / Chinese Astronomy and Astrophysics 30 (2006) 89"101
integrals, the system is integrable, and the phase orbits are quasi-periodic ordered orbits on the so-called KAM torus, the intersection of the five-dimensional hypersurfaces defined by the independent isolating integrals[28,24]. This integrable system generally has only 3 independent isolating integrals, but when there are 2 or 3 commensurable characteristic frequencies, there should exist 4 or 5 independent isolating integrals. On the other hand, if the number of independent isolating integrals possessed by the system is less than 3, then the system is non-integrable, which means, after the dynamic parameter passes a certain critical value, most of the orbits are chaotic orbits that are seemingly randomly distributed in the phase space and are no longer confined to the intersection of the five-dimensional hypersurfaces defined by the independent isolating integrals. To show more clearly that the isolating integral is a constraint on the phase space, we now take the two-dimensional linear oscillator as an illustrative example. This model is 1 2 1 22 H = ~(p~ + p ~ ) + ~(w~q~ +w2q2),
(2)
where wl and w2 are two fixed frequencies. It can be regarded as the composition of two independent one-dimensional linear harmonic oscillators along the ql and q2 axes, 1
Hi=~(pi
2
22 +wiqi)
(i=1,2).
(3)
Actually, H, H1, a n d / / 2 are all isolating integrals and we can write out a further isolating integral,
=/-/2 - H , .
(4)
Evidently, of the four isolating integrals, only two are independent. This means that the system is integrable and the phase space is two-dimensional. This fact can be further clarified by a canonical transformation to angle-action variables. Let ¢1 and ¢2 be the angle variables and J1 and -12, the action variables. Then the system (2) is transformed into H = l(waJl
+ "~2J2) •
(5)
-/2 = constant,
(6)
022 992 = ~ * t q- ~ 2 ( 0 ) ,
(7)
Hence, we have J1 = constant, Odl 991 --~-~ * t -}- ~ 1 ( 0 ) ,
and hence the time evolutions of the coordinates and momenta, qs = cs
sin(27c99~)(s = 1, 2),
p~ = c s w ~
cos(27r99~),
(8) (9)
where the constant ca = ~ (s = 1, 2). The two invariants of motion J1 and J2 of Eq.(6) reflect the two independent one-dimensional harmonic oscillator energy integrals H 1 = w l J ~ / 2 7 c a n d / / 2 = w2J2/2~r, which constrain the phase orbit of the particle to lie on a two-dimensional torus made up of the two tori determined respectively by ¢1 and ¢2 and
LIU Fu-yao et al. / Chinese Astronomy and Astrophysics 30 (2006) 8~101
91
therefore J1 and J2 are two isolating integrals. Perhaps, someone might derive from the two equations of (7) another invariant of motion, ( ~ 1 -- ~ 1 ( 0 ) ) / ( ~ 2
-- ~ 2 ( 0 ) )
:
~dl//-d2,
(10)
and assert that there always exists another isolating integral independent of J1 and J2. This conclusion seems to be true. We now discuss this problem separately for two cases. (1) The two frequencies wl and w2 are commensurable, i.e. if wl/w2 is equal to a rational number n / m , (n and m relative primes). Without loss of generality, we can let wl = 1 / m and w2 = 1/n. Then, in conjunction with Eq.(10) we have 27~m~1 - 27rncf12
=
2~(m~1 (0) - neff2(0))
=
constant.
(11)
Eq.(11) is indeed another isolating integral independent of J1 and J2. Its specific expression is derived as follows: Write the solutions (8) and (9) in complex form,
P-£ + iqs = cse i ' 2 ~
(i = x/~-l; s = 1, 2).
(12)
~ds
From this we can write
( ~Pl + iql ) ' ~ = cTe ~'2~'~1 P2 + iq2
(13)
= c~e i2'~n~'2 .
(14)
After dividing the left and right sides of Eq.(13) respectively by the left and right sides of Eq.(14), we have
z = ( p l / w l + iql)m(p2/w2 + iq2) - n = cr~e2 nei'27r(m~1(O)-n~2(O)).
(15)
It is easily seen t h a t E q . ( l l ) is the argument of the complex number z, which is an invariant of motion. The invariant point z on the complex plane reflects the common part of the two tori defined by ¢1 and ¢2. T h a t is to say, the phase paths on the two-dimensional torus form a closed and non-dense set of periodic orbits. The real and imaginary parts of z, Rez and Imz, are each an isolating integral. However, there is the following algebraic relation among H , ~, Rez and Imz: (Rez) 2 + (Imz) 2
= m2mn-2n(H
-
~ ) m ( H ~- ~1) - n ,
(16)
therefore there are only three independent isolating integrals. The specific forms of Rez and Imz for some given values of m and n now follow. When wl = w2 = w, the system (2) is isotropic. For simplicity, let rn = n = 1 in Eq.(15). Then we have (17) Rez = PiP2 + w2 ql q2 , (18)
Imz = qlP2 - q2pl • A
A
A
Rez and Rez differ only by a constant factor and so do Imz and Imz. Imz is none other than the angular momentum integral that reflects the axial symmetry of the system. Write I A $1 = ~-~wRez,
(19)
92
LIU Fu-yao et al. / Chinese Astronomy and Astrophysics 30 (2006) 87-101
--
2c~
'
s3 =
(2o) (21)
Then $1, $2 and $3 satisfy the relation { S i , S j } = ~ i j k S k ( i , j , k = 1,2,3;~jk the LeviCivita tensor)[ 3]. Thus, $1, $2 and $3 generate the symmetry group SO(3) of rotation transformation in the three-dimensional space. When w, ~ w2, the system (2) is anisotropic. Take m = 2 and n = 1 as an example. We obtain from Eq.(15), p2(4p~ - q~) + 4qlq2P2,
(22)
• ~ = Im"----z---- q2(q 2 - 4p~) + 4 p l P 2 q l .
(23)
~IJl = Re-'~ =
Since the system has no symmetry, neither ¢1 nor ¢2 can be the angular momentum integral, and ¢, ¢1 and ~2 can no longer generate the symmetry group SO(3). (2) When Wl/W2 is an irrational number, though Eq.(10) always exists, it can not be transformed into the form of Eq.(11), so the phase orbits on the two-dimensional torus are everywhere dense, indefinitely circulating, non-closed, quasi-periodic orbits. T h a t is to say, Eq.(10) imposes no constraints at all on the region of motion. So, Eq.(10) can not be regarded as an isolating integral: it is a non-isolating integral. Thus, from the mathematical point of view the distinction between isolating and nonisolating integrals is not obvious. From the physical point of view, however, a non-isolating integral is meaningless while an isolating integral imposes a constraint in the phase space. The commonly called integrals are isolating integrals. It needs be pointed out that the third isolating integral whose existence in system (2) was the subject of our inquiry in the foregoing cases (1) and (2), is completely different from the third isolating integral sought by H6non et al. [2a] in their study of qualitative dynamics. Whether or not the system (2) possesses a third isolating integral, it always has ordered orbits, but the distribution of the phase orbits is somewhat different. If the system (2) is isotropic, then the orbit in the configuration space is an ellipse with center at the coordinate origin. Therefore, so long as the semi-major axis, eccentricity and the azimuth of the major axis are known, the ellipse is determined. Thus to describe the plane ellipse centered at the origin 3 parameters need be known, i.e., 3 isolating integrals are required. If only the semi-major axis and eccentricity are known (equivalent to two isolating integrals) and the azimuth of the major axis is not fixed (equivalent to the non-existence of a third isolating integral), then the orbits will understandably be densely distributed. Thus, the existence of this kind of third isolating integral is also of significance. Namely, whether it exists or not will decide whether the phase orbit is closed and dense or not and whether it is periodic or quasi-periodic. On the other hand, in the context of H6non and Heilies, whether a third isolating integral exists or not determines whether or not the system is integrable, and hence whether the orbits will be ordered or chaotic. The axially symmetric autonomous Hamiltonian system with 3 degrees of freedom possesses the energy integral and the angular momentum integral and if there exists a third isolating integral, then the system is integrable; if not, it is non-integrable and it will possibly generate chaos. Thus, the existence of the third isolating integral is
LIU Fu-yao et al. / Chinese Astronomy and Astrophysics 30 (2006) 87-101
93
the indispensable condition of the integrability of the system. For the autonomous H~nonHeilies system with two degrees of freedom, there exists the energy integral but not the angular momentum integral, so the search for a third isolating integral becomes a search for a second isolating integral in this case. Under the condition that the H~non-Heilies system is so weak as to be non-integrable, it is difficult to search for a second isolating integral by analytical means. H~non and Heilies attempted to observe whether or not the KAM torus was destroyed in the Poincar~ section, thereby inferring whether or not there exists a second isolating integral. Now, we would like to lay stress on the following point: in an integrable, autonomous Hamiltonian system with n degrees of freedom, there exist at least n independent isolating integrals, but only n of these are mutually involutive. For example, the total energy integral H in Eq.(2) and S1, $2 and $3 in Eqs.(19)--(21) are all involutory but no two among the integrals $1, $2 and $3 are involutory. The meaning of isolating integral when discussing integrability in some books and documents may be different from that given in References [23, 24]. The isolating integrals referred to by Lichtenberg et al. [2] when they pointed out that "the essential condition of the integrability of the Hamiltonian system with n degrees of freedom is the existence of n independent isolating integrals" , are actually mutually involutory. From another direction, integrability of a system can also be understood in terms of the Lyapunov exponent. It is known that the autonomous Hamiltonian system with n degrees of freedom possesses 2n Lyapunov exponents of the form -t-Ai(i = 1 , . . . n) and Casati et al. proved that every isolating integral corresponds to two zero Lyapunov exponents [25]. Thus, for the integrable system with n independent isolating integrals, its Lyapunov exponents are all zero, and its orbit is generally quasi-periodic, which is the basis for the absence of chaos in integrable systems. Here, the n independent isolating integrals must still satisfy the mutually involutive relation. We now proceed to investigate the isolating integrals in the non-linear, isotropic Kepler two-body system. 2.2 I s o l a t i n g i n t e g r a l s in t h e K e p l e r t w o - b o d y m o t i o n The Hamiltonian for the Kepler two-body motion is 1
p
H = 2P •P-r' -
(24)
where r is the distance between the two bodies. H itself is an integral, called the energy integral. Because of spatial isotropy, there is rotational symmetry, hence conservation of angular momentum, i.e. L = r x p, (25) and simultaneously there exists the Runge-Lenz vector [26] R = p x L-
r
#-. (26) r The existence of these three integrals means the Kepler two-body system is integrable. For the Kepler two-body motion in a plane, the Hamiltonian has the form g = (p~ + p v2) / 2 - , / V / ~
+ y2 ,
(27)
and Eqs.(25) and (26) take the specific forms, L = (0, O, xpy - YPz),
(28)
94
LIU Fu-yao et al. / Chinese Astronomy and Astrophysics 30 (2006) 8~101
R =
p y(Xpy - Ypx)
~ , p#x x(yp~
-#Y -
- Xpy)
, 0) .
(29)
The nonzero components in Eqs.(28) and (29) are, in order, written as L z , R~ and R u. Now carry out the Levi-Civita canonical transformation
X ----2Ultt2 ,
1 u~p2 + u2Pl
Pz-- 2
u~+u~
'
(30)
1 U l p l - u2p2
Pv= 2
(31)
u~+u~
Then Eq.(27) is equivalent to y _ p~ + p~ - 8~
8(~1~ + u~)
(32)
Let dt = 4(Ul2 + u22)d~-, and the new Hamiltonian system obtained by a canonical extension of Eq.(32) is r
= 4(~I + ~ ) ( ~
+ po)
= 1(p12 +p~) + 4po(u~ + u 2) - 4#, 2
(33)
where P0 = - E (E is the energy of the system (27) and for the elliptical motion E < 0) and F - 0. F stands for the isotropic linear oscillator with frequency w = 8v/8~. This is also called the normalization method. Via this kind of processing, the Lyapunov unstable Kepler orbit is changed into the stable linear isotropic harmonic oscillator orbit. At the same time, Lz, R~ and R y formally become
1
(34)
L z = ~(u2pa - u l p 2 ) ,
1
Rx = - ~ (PIP2 + 8poulu2),
1 [ ( 12
G = ~
12
5p2 + @ 0 ~ ) - (5pl + @ o ~
(35)
)]
•
(36)
These three isolating integrals are extremely similar to Eqs.(19)-(21): we only have to multiply them with suitable constant factors to obtain the same forms, £,z = - L z ,
(37)
~ ,
(38)
2G
(39)
~_
95
LIU Fu-yao et al. / Chinese Astronomy and Astrophysics 30 (2006) 8~101
Thus, we can transform the 3 isolating integrals Lz, Rx and Ry of the system (27) into the isolating integrals Lz, Rx and Ry of the isotropic linear oscillator (33), and generate the symmetry group SO(3) of rotational transformation in three-dimensional space of the system (33). We now seek a numerical method that will conserve this symmetry group.
3. M A I N T A I N I N G T H E C O N S E R V E D Q U A N T I T I E S OF T H E K E P L E R MOTION WITH THE TRAPEZOIDAL FORMULA
For the first order ordinary differential equation
= f ( t , x),
(40)
one may resort to the linear k-step implicit Newton-Coates quadrature formula, k ~gn+k : 2gn 3- h E
Cikf n+i'
(41)
i=0
where
(_l)k_i ~ok k
Cik -- iV~(k=~) !
I](s - j)ds, j¢i
fn+i = f(tn+i, xn+i),
to obtain numerical solutions, with h the integration step and tn+i = tn 3- ih. Here, a few Newton-Coates integrators are enumerated: h
X n + l -- Xn : - 2 ( f n + f n + l ) ,
k : 1,
]¢ ----2,
k = 3, k=4,
~lgn+3
~n+2 -- ggn ---- h ( f n
--
3- 4 f n + l 3- f n + 2 ) ,
Xn = 3--hs(fn + 3fn+ 1 + 3 f n + 2 3- f n + 3 ) ,
Xn+4 --X~ = -~2h(T f n 3- 3 2 f n + l 3- 12fn+2 3- 3 2 f n + 3 3- 7 fnJc4 ) .
(42)
(43)
(44) (45)
They are called, in order, the second order single-step trapezoidal formula, the fourth order two-step Simpson formula, the fourth order three-step Newton formula and the sixth order four-step Cotes formula. According to Feng Kangls theory [15], it is not difficult to judge that the linear kstep implicit Newton-Coates quadrature formula (41) is symplectic. Feng Kang/s theory is described as follows: For the linear k-step algorithm, k
k
E O~iXn+i = h E ~ifn+i i=O
i=O
(46)
96
LIU Fu-yao et al. / Chinese Astronomy and Astrophysics 30 (2006) 8~101
used to solve linear Hamiltonian systems to be symplectic, it is necessary and sufficient that the two characteristic polynomials of Eq.(46), k
p(~) = ~-~' ai~ ~,
k
a(~) = ~-~.~
i=0
(47)
i=0
should satisfy the relation ¢(~) = - ¢ ( 1 / ~ ) (~(~) = p(~)/a(~)), or that p(~) is antisymmetric while a(~) is symmetric. Therefore we know that Eq.(41) is a linear multi-step symplectic algorithm. What is worth noting is that the symplectic property shown by the linear multi-step formula (41) is tenable only when the method is numerically stable and when the solutions of the difference equation of the linear multi-step method are determined solely by the principal root, free of any effect from the spurious roots [16,17]. Excepting the trapezoidal formula (42) which has a comparatively large and relatively stable region, all other NewtonCoates integration formulae have zero stability, i.e. the stable region is only limited to certain segments symmetrical about the origin on the imaginary axis of the complex plane. For linear systems, so long as the step length is proper, the stability of these linear multistep methods is guaranteed. When they are applied to non-linear systems, rather strong numerical instability will emerge. 3.1 C o n s e r v a t i o n of t h e s y m m e t r y g r o u p of t h e i s o t r o p i c t w o - d i m e n s i o n a l barmonic oscillator The second order single-step trapezoidal formula (42) has a much better stability than other Newton-Coates integration methods and since it is also a sing-step method, it has obvious practical values. But in practice, it has not been given much attention, because it is an implicit method, which generally requires iteration, and so has a lower efficiency than an explicit algorithm. However, when applied to a linear systems an implicit method can be converted into an explicit one. First, for a brief presentation of the problem we shall analyze the harmonic oscillator formula (3) (omitting the subscript i). The differential equation of motion here is (~ = p,
~5= _w2q.
(48)
Its analytic solution is q(t) = qo cos(wt) + ( 1 / w ) p o sin(wt), p(t)
= - w q o sin(wt) + P0 cos(wt).
(49)
Note that Eqs.(49) are only formally different from Eqs. (8) and (9). Applying the trapezoidal formula (42) on Eq.(48), we obtain numerical solutions of the form, qn+l = qn cosc~ + (1/W)pn sin(~, P~+I = -Wqn sin (~ + pn cos c~. Here cosa = ( 4 - w2h2)/(4 +w2h2),
(50)
LIU Fu-yao et al. / Chinese Astronomy and Astrophysics 30 (2006) 87-101
97
sinc~ = 4wh/(4 +w2h2),
~ hf~(h) f~(h) = 4w/(4 - w2h 2) ~-, w + w3h2/4
(51)
It is easily verified that dqn+l A d p n + l = dqn A d p ~ , i.e. Eq.(50) is area-preserving or is a symplectic mapping t. The numerical solution Eq.(50) and the analytic solution Eq.(49) are completely identical in form, but the m a p p e d frequency ft(h) is evidently not equal to the proper frequency w of the system. Therefore, the numerical solution will manifest a linearly increasing along-track error but will not introduce any artificial dissipation. It is not difficult to get H(qn+~,pn+~) - H(q~,pn), showing that the trapezoidal formula strictly conserves the energy of the harmonic oscillator, or produces no truncation error in the energy. This is different from the conservation of energy by the usual symplectic algorithm, which always contains the truncation error in energy which, however, always fluctuates about the original energy, thereby maintaining the energy in a sense, and it m a y be said that the ordinary symplectic algorithm maintains an energy integral in an altered form. Similarly, let us solve numerically the isotropic two-dimensional harmonic oscillator (33) (frequency w = 8 v ~ ) by roans of the trapezoidal formula (42). Then from the result obtained we easily deduce F - 0, and t h a t the 3 isolating integrals Lz, Rx and /~y in Eqs.(37)--(39) are strictly conserved, i.e., the s y m m e t r y group SO(3) they generate is strictly conserved. Then, the isolating integrals in Eqs.(27)--(29) are strictly conserved by the trapezoidal formula (42). 3.2
Conservation of the 5 Kepler orbital elements For the Kepler two-body motion in three-dimensional space, we start with a total of 7 isolating integrals i.e., 1 energy integral, 3 components of the angular m o m e n t u m integral and 3 components of the Runge-Lenz vector (Eqs.(24)-(26)). But since R . L = 0, R and the radial unit v e c t o r / ~ of the periastron point satisfies the relation z
r
--
1 (p.p)r_(r.p)p_#r e#
~ -e
--
=R/(e#)'
--
e#
r - ~ -e
--
--
e#
p= (52)
the Runge-Lenz vector is parallel to the periastron vector and its modulus satisfies the relation IRI = e , (e = IRI/p). (53) Hence there are only 5 independent isolating integrals. (If we use the canonical Delaunay variables, then there will be 5 cyclic coordinates in the new Hamiltonian, producing 5 cyclic integrals and making the existence of 5 independent isolating integrals much more obvious). It is rather obvious t h a t to specify an ellipse centered at the origin in three-dimensional space we need 5 independent parameters, namely, the semi-major radius a, eccentricity e, orbital inclination i, longitude of node f~ and periastron angle w, and hence 5 independent isolating integrals. 4 of the 5 orbital elements other t h a n e can be expressed in terms of the t Note that there exists still the truncation error in the numerical solution which shows up in the mapped frequency ft(h).
98
LIU Fu-yao et al. / Chinese Astronomy and Astrophysics 30 (2006) 87-101
integrals (24)-(26) as follows[26]: a = -#/(2H),
(54)
i = arccos R z ,
(55)
12 = arctan(- Rx / Ry ) ,
(56)
w = arcsin(/3z/sin i) = arcsin[Rz/(e# sin i)].
(57)
Of course, to place ~ and w in their right quadrants, we need the signs of certain trigonometric functions. Minesaki et al. [22] adopted the KS transformation to transform the non-linear Kepler two-body system into the four-dimensional isotropic harmonic oscillator, on which they then applied the trapezoidal formula; they then recovered the original system by inversion, all the while strictly conserving the energy integral, angular momentum integral and Runge-Lenz vector. From this we can conclude that the trapezoidal formula can strictly conserve the 5 Kepler orbital elements a, e, i, ~ and w in the light of Eqs.(53)-(57). Here, we should emphasize that the 5 orbital elements only specify the ellipse in relation to the coordinate origin in three-dimensional space; to specify a particular point on the ellipse, we need a further parameter, the mean anomaly M, say. Now, let us see the effect of the error in mean anomaly M caused by the trapezoidal formula. In Reference [22] a fixed step size is adopted for the physical time t; it is changed to a variable step size for the new coordinate time after the KS transformation. Since a is strictly conserved by the trapezoidal formula, the mean motion n is also strictly conserved. For the physical time covered in one step, the increment of M is A M ----n. At, showing that no truncation error on M is produced by the trapezoidal formula. In contrast, for the usual symplectic algorithm, because there exists a truncation error A n in n, we have ~ = n + A n , and the increment of M is A M = ~ . A t = n . A t + A n . A t , where An. At is the local truncation error in M produced by the usual symplectic algorithm. This shows that the error in M caused by the usual symplectic algorithm induces a longterm variation. This is another aspect where the trapezoidal formula is superior to the usual symplectic algorithm.
4. C O U P L E D
HARMONIC
OSCILLATOR
For the anisotropic linear oscillator, although the original system has the proper frequencies Wl/W2 = n / m , the mapped frequencies (Eq.(51)) obtained through the trapezoidal formula have ~tl(h)/~2(h) ~ n / m . Thus, the trapezoidal formula does not conserve other isolating integrals than the energy integral (or some integral associated with it), that is, the isolating integrals (22) and (23)). Let us take up the coupled harmonic oscillator once again. The coupled harmonic oscillator is H = ~ l ( p ~ +p22 ) + 2lw2r o(.ql2 +q22) - 6qlq2
(58)
with ~ the coupling parameter. Make the coordinate transformation
ql-
q2-
(59)
LIU Fu-yao et al. / Chinese Astronomy and Astrophysics 30 (2006) 87-101
99
and obtain the new Hamiltonian system
1 2 +p~) + ! {022~2 +w~/2) 2 ~ 1~ ,
H = ~(p~
(60)
021 ~- W0v / f l - £/0.)02 ,
(61)
w2 = wo V/I + e/w~ .
(62)
Eq.(60) stands for the anisotropic linear oscillator. When c/w~ << 1, Wl ~ w0[1-e/(2w02)] and ~2 ~ w0[1 + c/(2w0z)]. This indicates t h a t no m a t t e r how small e/w~ is, within an accuracy greater t h a n the accuracy of the computer, the trapezoidal formula no longer conserves the isolating integrals of the coupled harmonic oscillator except the energy integral or any integral associated with it. Thus, the condition of strict conservation of the relevant isolating integrals by the trapezoidal formula is very stringent. Even if the system is a linear system, as soon as it is anisotropic, then the trapezoidal formula ceases to conserve strictly all the isolating integrals except the energy integral or any integral related to it. The Kepler two-body system excepted, all other nonlinear integrable autonomous Hamiltonian systems with n degrees of freedom (such as the one Consisting of a central gravitational field and an additional uniform force field [1] or the problem of two fixed centers) can not, in general, be transformed into an isotropic harmonic oscillator. If there exist more t h a n n independent isolating integrals, the trapezoidal formula generally can strictly conserve at most n such integrals.
5. C O N C L U D I N G
REMARKS
The linear autonomous Hamiltonian system with n degrees of freedom certainly has n mutually involutive independent isolating integrals, and these independent isolating integrals are strictly conserved by the trapezoidal formula. If there exists commensurability among some or all of the characteristic frequencies, then the number of independent isolating integrals can exceed n, but can not exceed 2n - 1. When the system is also isotropic, it possesses 2 n - 1 independent isolating integrals, and all of t h e m can be strictly conserved by the trapezoidal formula. If all the characteristic frequencies of the system are commensurate and the system is isotropic then even though there exist 2n - 1 independent isolating integrals, none of these are conserved by the trapezoidal formula except the energy integral (or some other integral related to it). The two-dimensional Kepler two-body system pertains to an isotropic nonlinear system and has the energy integral, the angular m o m e n t u m integral and the Runge-Lenz vector. In the case of plane motion, among these three quantities there are a total of 4 nonzero isolating integrals but only 3 of them are independent ones. These 3 independent isolating integrals make up the s y m m e t r y group SO(3) of the system, and via the Levi-Civita transformation it cam be transformed into the s y m m e t r y group of the two-dimensional isotropic harmonic oscillator, which is strictly conserved by the trapezoidal formula. In the case of motion in three-dimensional space, the three n a m e d quantities have in total 7 isolating integrals of which only 5 are independent, as the Runge-Lenz vector is parallel to the the periastron
100
LIU Fu-yao et al. / Chinese Astronomy and Astrophysics 30 (2006) 8"7101
point vector. Through the KS transformation the system can be transformed into the four-dimensional isotropic harmonic oscillator (in topological language, the two- and threedimensional Kepler two-body systems are respectively homeomorphic to the two- and fourdimensional isotropic harmonic oscillators). The trapezoidal formula can likewise strictly conserve the three named quantities (the energy and angular momentum integrals and the Runge-Lenz vector), and this fact is reflected, in particular, in the 5 Kepler orbital elements a, e, i, gt and w. As for other nonlinear integrable autonomous Hamiltonian systems, generally speaking, they can not be transformed into isotropic harmonic oscillators, and the trapezoidal formula can not in general conserve isolating integrals with forms similar to the form of the Runge-Lenz vector. It needs to be emphatically pointed out that although it is seen analytically that the trapezoidal formula strictly conserves the energy integral, angular momentum integral and Runge-Lenz vector as well as the 5 Kepler orbital elements a, e, i, ~t and w, numerical results obtained by the computer always have errors that increase linearly with the integration time. This is because freedom from the truncation error makes the rounding error the principal source of error. In passing, it may be mentioned that the equations of stellar orbital motion and of the relativistic perihelion precession can be simply and conveniently derived by taking advantage of the Runge-Lenz vector [27]. ACKNOWLEDGEMENT The second author of this article discussed the concepts of isolating and non-isolating integrals with Prof. H U A N G Tian-yi of Department of Astronomy of Nanjing University and derived great benefits from the discussion. We wish to express our indebtedness to Prof. Huang. References 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18
Arnold V.I., Mathematical Methods of Classical Mechanics (translated by Vogtmann K. and Weinstein A.). New York: Springer-Verlag, 1978, 100 LichtenbergS.J. Lieberman M.A., Regular and Chaotic Dynamics. Springer-Verlag, 1983, 1 GoldsteinH., Classical Mechanics. Addison-Wesley Publishing Company, 1950, 55-300 Wu X., Huang T.Y, Wan X.S., ChASzA, 2003, 27(1). 114 Wu X., Huang T.Y., Wan X.S., AcASn, 2002, 43(4), 391 Liu F.Y., Wu X., Lu B.K., ChA&A, 2005, 29(1), 92 Liu F.Y., Wu X., Lu B.K., AcASn, 2004, 45(4), 402 Wu X., Huang T.Y., Wan X.S., ChA&A, 2005, 29(1), 81 Wu X., Huang T.Y., Wan X.S., AcASn, 2004, 45(3), 310 WisdomJ., Holman M., AJ, 1991, 102, 1528 YoshidaH., CeMDA, 1993, 56, 27 AscherU.M., Numerical Algorithms, 1997, 14, 1 YoshidaH., Phys. Lett. A, 2001, 282, 276 YoshidaH., CeMDA, 2002, 83, 355 Feng K., Proceedings of the Annual Meeting on Computational Mathematics in Tianjin. Tianjin, Nankai University Press, 1990, 1 Huang T.Y., Wang C.B., Zhao Z.Y., ChA&:A, 1998, 22, 101 Huang T.Y., Wang C.B., Zhao Z.Y., AcASn, 1997, 38, 278 Chiou J.C., Wu S.D., J. Chem. Phys., 1997, 107(17), 6894
LIU Fu-yao et al. / Chinese Astronomy and Astrophysics 30 (2006) 87-101
19 20 21 22 23 24 25 26 27
101
Zhu W., Zhao X., Tang Y., J. Chem. Phys., 1996, 104(6), 2275 Kalogiratou Z., Simos T.E., J. Comput. ~ App. Math., 2003, 158, 75 Minesaki Y., Nakamura Y., Phys. Lett. A., 2002,306, 127 Minesaki Y., Nakamura Y., Phys. Lett. A., 2004, 324, 282 Hdnon M., Heiles C., AJ, 1964, 69, 73 http://wwwgoogle.com. Annu. Rev. Astron. Astrophys., 1982, 20, 399: Casati G., Chirikov B.V., Ford J., Phys. Lett. A., 1980, 77, 91 Brouwer D., Clemence G.M., Methods of Celestial Mechanics, New York: Academic Press, 1961, 15 Weinberg S., Gravitation and Cosmology, New York: John Wiley, 1972, 165
ELSEVIER Chinese Astronomy and Astrophysics 30 (2006) 87 101
A Numerical Method Runge-Lenz
that Conserves Vector t *
the
LIU Fu-yao 1,4 W U Xin 2'3 LU Ben-kui 1'4 1purple Mountain Observatory, Chinese Academy of Sciences, Nanjing 210008 2School of Science, Nanchang University, Nanchang 330047 3Department of Astronomy, Nanjing University, Nanjing 210093 4National Astronomical Observatories, Chinese Academy of Sciences, Beijing 100012
A b s t r a c t An exhaustive discussion is carried out on isolating integrals and the trapezoidal formula which can conserve the Runge-Lenz vector. An isolating integral is an invariant that restricts the region of particle motion. The autonomous integrable Hamiltonian system with n degrees of freedom has only n mutually involutive independent isolating integrals, and the existence of other isolating integrals is meaningful to the particle motion. In the Kepler two-body system there exist the energy integral, the angular momentum integral and the RungeLenz vector. These correspond to 3 independent isolating integrals in the case of plane motion, and to 5 in the case of space motion. In the former, the integrals makes up the symmetry group SO (3) of the system, which can be transformed into the symmetry group of the two-dimensional isotropic harmonic oscillator through the Levi-Civita transformation, which is accurately conserved by the trapezoidal formula. On the other hand, in the case of space motion, the strict conservation of the energy and angular momentum inegrals and the Runge-Lenz vector by the trapezoidal formula is manifested in the 5 Kepler orbital elements a, e, i, ft and w. K e y w o r d s : celestial mechanics--methods: numerical
t Supported by National Natural Science Foundation Received 2004 11 01; revised version 2004-11-22 * A translation of Acta Astron. Sin. Vol. 46, No. 3, pp. 294-306, 2005 0275-1062/06/S-see front matter © 2006 Elsevier B. V.All rights reserved. DOI: 10.1016/j.chinastron.2006.01.007
88
LIU Fu-yao et al. / Chinese Astronomy and Astrophysics 30 (2006) 87-101
1. I N T R O D U C T I O N In terms of linear and nonlinear relations in kinetic differential equations, Hamiltonian systems can be divided into two kinds, linear and nonlinear systems. The linear systems and some of the nonlinear systems belong to integrable systems, possessing the dynamic properties of ordered motion free of chaotic behaviour. The necessary and sufficient condition for an autonomous Hamiltonian system with n degrees of freedom to be integrable is that there exist n mutually involutive, independent isolating integrals [1,2]. Of course, some integrable Hamiltonian systems of n degrees of freedom possess more than n independent isolating integrals. For a linear system of n degrees of freedom with frequency commensuration, there exist 2n - 1 independent isolating integrals;-- one such example is the isotropic autonomous harmonic oscillator of n degrees of freedom. For the Kepler two-body system which pertains to an isotropic system, there exist the energy integral, the angular momentum integral and the Runge-Lenz vector. The first two respectively reflect invariance with respect to time translation and space rotation. In the case of plane motion, and all three represent independent isolating integrals, while in the case of space motion there are only 5 independent integrals [a] out of a total of 7 integrals (i.e., 1 of energy, 3 of the components of angular momentum and 3 of the components of the Runge-Lenz vector). A good numerical method should have a good computational efficiency, a high numerical accuracy as well as a good numerical stability and, especially, can conserve the geometrical topological structure and proper invariant manifold of the original system. Examples are the symplectic algorithm[ 4-7], the stabilization method and the method of reduction of order of equation Is,9]. Unquestionably, the often used symplectic algorithms such as the first order symplectic Euler method and the second order Wisdom-Holman leap-frog symplectic mapping, etc., are able to maintain the energy integral[ 1°], but this does not mean that the symplectic algorithm strictly conserves the energy integral. As a m a t t e r of fact, there always exists the truncation error In] when the Hamiltonian system obtained from the numerical results is compared with the original Hamiltonian system. However, the error in the Hamiltonian function has no long-term variation, so that the integration always proceeds near the plane of energy. Thus, the symplectic algorithm conserves an altered form of the energy integral, that is to say, the symplectic algorithm is a method which automatically satisfies the energy constraint condition. From this angle it is seen that the symplectic algorithm generally needs no stabilization treatment [12]. Besides, for the isotropic harmonic oscillator and the Kepler two-body system, Yoshida found that the symplectic algorithm can conserve the angular momentum integral in an altered form [13,14]. However, the symplectic algorithm, generally speaking, can not conserve all the isolating integrals. Yoshida pointed out that when the symplectic algorithm is applied to a linear anisotropic system with frequency commensuration, no other integrals in their altered forms exist except the energy integral in its altered form [131. He made further discussions and pointed out that the symplectic algorithm is unable to conserve the proper Runge-Lenz vector of the Kepler two-body system [14]. However, there is another sort of symplectic algorithms which can strictly conserve the energy integral and other integrals of motion. The simplest trapezoidal formula in the Newton-Coates family of integration as a linear multi-step implicit symplectic algorithm can accurately conserve the energy of the linear oscillatorI15-171. Here, it should be emphasized that the linear multi-step symplectic algorithm and the symmetrized algo-
LIU Fu-yao et al. / Chinese Astronomy and Astrophysics 30 (2006) 87-101
89
rithm are actually approximate symplectic structure preserving algorithms, i.e. subject to the conditions that the algorithm is stable and that the characteristic equation is determined by the principal root while the effect of spurious roots can be neglected [15-17]. For the linear Hamiltonian system, a sort of multilayer symplectic integrator in the shape of the trapezoidal formula can be constructed by making use of the Newton-Coates formula and the implicit midpoint formula [ls-2°]. As for the nonlinear system, Minesaki and Nakamura studied the feasibility of the application of the trapezoidal formula. Using the Levi-Civita or KS transformation, they transformed the nonlinear Kepler two-body system into the isotropic harmonic oscillator, to which they applied the trapezoidal formula. They then worked back to the original system, thus with the energy integral, the angular momentum integral and the Runge-Lenz vector are strictly conserved [21,221. Considering that the deduction by Minesaki and Nakamura of strict conservation of the Runge-Lenz vector with the trapezoidal formula is somewhat tedious, we seek in this paper a simpler derivation following a detailed exploration. The integrals of the Kepler two-body plane motion make up the symmetry group SO (3) of the system, which can be transformed into the symmetry group of the two-dimensional isotropic harmonic oscillator via the Levi-Civita transformation and the latter symmetric group can be strictly conserved by the trapezoidal formula. This way of derivation is comparatively simple. At the same time we shall make an in-depth exploration of isolating integral and the condition of integrability. We shall also examine the numerical accuracy of the six Kepler orbital elements. Five of the elements, i.e. a, e, i, ~ and ~, should be constants and we shall investigate whether they are strictly conserved by the trapezoidal formula. We shall also explore the condition under which the isolating integrals are conserved by numerical integration of the nonlinear integrable system including the two-dimensional coupled harmonic oscillator.
2. S Y M M E T R Y G R O U P OF T H E K E P L E R T W O - B O D Y M O T I O N For later convenience, we state in this section the concepts of isolating and non-isolating integrals as well as the symmetry group of the Kepler two-body plane motion, which can be transformed into that of the two-dimensional isotropic harmonic oscillator.
2.1 Isolating and non-isolating integrals Research on dynamical systems generally revolves around such concepts as isolating integral, non-isolating integral, integrability and non-integrability, etc. An isolating integral is the very invariant that restricts the region of motion of the particle. For an autonomous Hamiltonian system H = H ( q , p), if a function (I) = (I)(q, p) (of the generalized coordinate and momentum) is such that its Poisson bracket with H satisfies the involutory relation,
{e, H} = 0,
(1)
then (P is said to be an isolating integral of H. An autonomous Hamiltonian system of 3 degrees of freedom has at least one isolating integral because H itself is an invariant of motion (called the energy integral). It has at most 5 independent isolating integrals, because the number of independent variables in the phase space is five, and each isolating integral is a five-dimensional hypersurface. When the system possesses 3 or more independent isolating
90
LIU Fu-yao et al. / Chinese Astronomy and Astrophysics 30 (2006) 89"101
integrals, the system is integrable, and the phase orbits are quasi-periodic ordered orbits on the so-called KAM torus, the intersection of the five-dimensional hypersurfaces defined by the independent isolating integrals[28,24]. This integrable system generally has only 3 independent isolating integrals, but when there are 2 or 3 commensurable characteristic frequencies, there should exist 4 or 5 independent isolating integrals. On the other hand, if the number of independent isolating integrals possessed by the system is less than 3, then the system is non-integrable, which means, after the dynamic parameter passes a certain critical value, most of the orbits are chaotic orbits that are seemingly randomly distributed in the phase space and are no longer confined to the intersection of the five-dimensional hypersurfaces defined by the independent isolating integrals. To show more clearly that the isolating integral is a constraint on the phase space, we now take the two-dimensional linear oscillator as an illustrative example. This model is 1 2 1 22 H = ~(p~ + p ~ ) + ~(w~q~ +w2q2),
(2)
where wl and w2 are two fixed frequencies. It can be regarded as the composition of two independent one-dimensional linear harmonic oscillators along the ql and q2 axes, 1
Hi=~(pi
2
22 +wiqi)
(i=1,2).
(3)
Actually, H, H1, a n d / / 2 are all isolating integrals and we can write out a further isolating integral,
=/-/2 - H , .
(4)
Evidently, of the four isolating integrals, only two are independent. This means that the system is integrable and the phase space is two-dimensional. This fact can be further clarified by a canonical transformation to angle-action variables. Let ¢1 and ¢2 be the angle variables and J1 and -12, the action variables. Then the system (2) is transformed into H = l(waJl
+ "~2J2) •
(5)
-/2 = constant,
(6)
022 992 = ~ * t q- ~ 2 ( 0 ) ,
(7)
Hence, we have J1 = constant, Odl 991 --~-~ * t -}- ~ 1 ( 0 ) ,
and hence the time evolutions of the coordinates and momenta, qs = cs
sin(27c99~)(s = 1, 2),
p~ = c s w ~
cos(27r99~),
(8) (9)
where the constant ca = ~ (s = 1, 2). The two invariants of motion J1 and J2 of Eq.(6) reflect the two independent one-dimensional harmonic oscillator energy integrals H 1 = w l J ~ / 2 7 c a n d / / 2 = w2J2/2~r, which constrain the phase orbit of the particle to lie on a two-dimensional torus made up of the two tori determined respectively by ¢1 and ¢2 and
LIU Fu-yao et al. / Chinese Astronomy and Astrophysics 30 (2006) 8~101
91
therefore J1 and J2 are two isolating integrals. Perhaps, someone might derive from the two equations of (7) another invariant of motion, ( ~ 1 -- ~ 1 ( 0 ) ) / ( ~ 2
-- ~ 2 ( 0 ) )
:
~dl//-d2,
(10)
and assert that there always exists another isolating integral independent of J1 and J2. This conclusion seems to be true. We now discuss this problem separately for two cases. (1) The two frequencies wl and w2 are commensurable, i.e. if wl/w2 is equal to a rational number n / m , (n and m relative primes). Without loss of generality, we can let wl = 1 / m and w2 = 1/n. Then, in conjunction with Eq.(10) we have 27~m~1 - 27rncf12
=
2~(m~1 (0) - neff2(0))
=
constant.
(11)
Eq.(11) is indeed another isolating integral independent of J1 and J2. Its specific expression is derived as follows: Write the solutions (8) and (9) in complex form,
P-£ + iqs = cse i ' 2 ~
(i = x/~-l; s = 1, 2).
(12)
~ds
From this we can write
( ~Pl + iql ) ' ~ = cTe ~'2~'~1 P2 + iq2
(13)
= c~e i2'~n~'2 .
(14)
After dividing the left and right sides of Eq.(13) respectively by the left and right sides of Eq.(14), we have
z = ( p l / w l + iql)m(p2/w2 + iq2) - n = cr~e2 nei'27r(m~1(O)-n~2(O)).
(15)
It is easily seen t h a t E q . ( l l ) is the argument of the complex number z, which is an invariant of motion. The invariant point z on the complex plane reflects the common part of the two tori defined by ¢1 and ¢2. T h a t is to say, the phase paths on the two-dimensional torus form a closed and non-dense set of periodic orbits. The real and imaginary parts of z, Rez and Imz, are each an isolating integral. However, there is the following algebraic relation among H , ~, Rez and Imz: (Rez) 2 + (Imz) 2
= m2mn-2n(H
-
~ ) m ( H ~- ~1) - n ,
(16)
therefore there are only three independent isolating integrals. The specific forms of Rez and Imz for some given values of m and n now follow. When wl = w2 = w, the system (2) is isotropic. For simplicity, let rn = n = 1 in Eq.(15). Then we have (17) Rez = PiP2 + w2 ql q2 , (18)
Imz = qlP2 - q2pl • A
A
A
Rez and Rez differ only by a constant factor and so do Imz and Imz. Imz is none other than the angular momentum integral that reflects the axial symmetry of the system. Write I A $1 = ~-~wRez,
(19)
92
LIU Fu-yao et al. / Chinese Astronomy and Astrophysics 30 (2006) 87-101
--
2c~
'
s3 =
(2o) (21)
Then $1, $2 and $3 satisfy the relation { S i , S j } = ~ i j k S k ( i , j , k = 1,2,3;~jk the LeviCivita tensor)[ 3]. Thus, $1, $2 and $3 generate the symmetry group SO(3) of rotation transformation in the three-dimensional space. When w, ~ w2, the system (2) is anisotropic. Take m = 2 and n = 1 as an example. We obtain from Eq.(15), p2(4p~ - q~) + 4qlq2P2,
(22)
• ~ = Im"----z---- q2(q 2 - 4p~) + 4 p l P 2 q l .
(23)
~IJl = Re-'~ =
Since the system has no symmetry, neither ¢1 nor ¢2 can be the angular momentum integral, and ¢, ¢1 and ~2 can no longer generate the symmetry group SO(3). (2) When Wl/W2 is an irrational number, though Eq.(10) always exists, it can not be transformed into the form of Eq.(11), so the phase orbits on the two-dimensional torus are everywhere dense, indefinitely circulating, non-closed, quasi-periodic orbits. T h a t is to say, Eq.(10) imposes no constraints at all on the region of motion. So, Eq.(10) can not be regarded as an isolating integral: it is a non-isolating integral. Thus, from the mathematical point of view the distinction between isolating and nonisolating integrals is not obvious. From the physical point of view, however, a non-isolating integral is meaningless while an isolating integral imposes a constraint in the phase space. The commonly called integrals are isolating integrals. It needs be pointed out that the third isolating integral whose existence in system (2) was the subject of our inquiry in the foregoing cases (1) and (2), is completely different from the third isolating integral sought by H6non et al. [2a] in their study of qualitative dynamics. Whether or not the system (2) possesses a third isolating integral, it always has ordered orbits, but the distribution of the phase orbits is somewhat different. If the system (2) is isotropic, then the orbit in the configuration space is an ellipse with center at the coordinate origin. Therefore, so long as the semi-major axis, eccentricity and the azimuth of the major axis are known, the ellipse is determined. Thus to describe the plane ellipse centered at the origin 3 parameters need be known, i.e., 3 isolating integrals are required. If only the semi-major axis and eccentricity are known (equivalent to two isolating integrals) and the azimuth of the major axis is not fixed (equivalent to the non-existence of a third isolating integral), then the orbits will understandably be densely distributed. Thus, the existence of this kind of third isolating integral is also of significance. Namely, whether it exists or not will decide whether the phase orbit is closed and dense or not and whether it is periodic or quasi-periodic. On the other hand, in the context of H6non and Heilies, whether a third isolating integral exists or not determines whether or not the system is integrable, and hence whether the orbits will be ordered or chaotic. The axially symmetric autonomous Hamiltonian system with 3 degrees of freedom possesses the energy integral and the angular momentum integral and if there exists a third isolating integral, then the system is integrable; if not, it is non-integrable and it will possibly generate chaos. Thus, the existence of the third isolating integral is
LIU Fu-yao et al. / Chinese Astronomy and Astrophysics 30 (2006) 87-101
93
the indispensable condition of the integrability of the system. For the autonomous H~nonHeilies system with two degrees of freedom, there exists the energy integral but not the angular momentum integral, so the search for a third isolating integral becomes a search for a second isolating integral in this case. Under the condition that the H~non-Heilies system is so weak as to be non-integrable, it is difficult to search for a second isolating integral by analytical means. H~non and Heilies attempted to observe whether or not the KAM torus was destroyed in the Poincar~ section, thereby inferring whether or not there exists a second isolating integral. Now, we would like to lay stress on the following point: in an integrable, autonomous Hamiltonian system with n degrees of freedom, there exist at least n independent isolating integrals, but only n of these are mutually involutive. For example, the total energy integral H in Eq.(2) and S1, $2 and $3 in Eqs.(19)--(21) are all involutory but no two among the integrals $1, $2 and $3 are involutory. The meaning of isolating integral when discussing integrability in some books and documents may be different from that given in References [23, 24]. The isolating integrals referred to by Lichtenberg et al. [2] when they pointed out that "the essential condition of the integrability of the Hamiltonian system with n degrees of freedom is the existence of n independent isolating integrals" , are actually mutually involutory. From another direction, integrability of a system can also be understood in terms of the Lyapunov exponent. It is known that the autonomous Hamiltonian system with n degrees of freedom possesses 2n Lyapunov exponents of the form -t-Ai(i = 1 , . . . n) and Casati et al. proved that every isolating integral corresponds to two zero Lyapunov exponents [25]. Thus, for the integrable system with n independent isolating integrals, its Lyapunov exponents are all zero, and its orbit is generally quasi-periodic, which is the basis for the absence of chaos in integrable systems. Here, the n independent isolating integrals must still satisfy the mutually involutive relation. We now proceed to investigate the isolating integrals in the non-linear, isotropic Kepler two-body system. 2.2 I s o l a t i n g i n t e g r a l s in t h e K e p l e r t w o - b o d y m o t i o n The Hamiltonian for the Kepler two-body motion is 1
p
H = 2P •P-r' -
(24)
where r is the distance between the two bodies. H itself is an integral, called the energy integral. Because of spatial isotropy, there is rotational symmetry, hence conservation of angular momentum, i.e. L = r x p, (25) and simultaneously there exists the Runge-Lenz vector [26] R = p x L-
r
#-. (26) r The existence of these three integrals means the Kepler two-body system is integrable. For the Kepler two-body motion in a plane, the Hamiltonian has the form g = (p~ + p v2) / 2 - , / V / ~
+ y2 ,
(27)
and Eqs.(25) and (26) take the specific forms, L = (0, O, xpy - YPz),
(28)
94
LIU Fu-yao et al. / Chinese Astronomy and Astrophysics 30 (2006) 8~101
R =
p y(Xpy - Ypx)
~ , p#x x(yp~
-#Y -
- Xpy)
, 0) .
(29)
The nonzero components in Eqs.(28) and (29) are, in order, written as L z , R~ and R u. Now carry out the Levi-Civita canonical transformation
X ----2Ultt2 ,
1 u~p2 + u2Pl
Pz-- 2
u~+u~
'
(30)
1 U l p l - u2p2
Pv= 2
(31)
u~+u~
Then Eq.(27) is equivalent to y _ p~ + p~ - 8~
8(~1~ + u~)
(32)
Let dt = 4(Ul2 + u22)d~-, and the new Hamiltonian system obtained by a canonical extension of Eq.(32) is r
= 4(~I + ~ ) ( ~
+ po)
= 1(p12 +p~) + 4po(u~ + u 2) - 4#, 2
(33)
where P0 = - E (E is the energy of the system (27) and for the elliptical motion E < 0) and F - 0. F stands for the isotropic linear oscillator with frequency w = 8v/8~. This is also called the normalization method. Via this kind of processing, the Lyapunov unstable Kepler orbit is changed into the stable linear isotropic harmonic oscillator orbit. At the same time, Lz, R~ and R y formally become
1
(34)
L z = ~(u2pa - u l p 2 ) ,
1
Rx = - ~ (PIP2 + 8poulu2),
1 [ ( 12
G = ~
12
5p2 + @ 0 ~ ) - (5pl + @ o ~
(35)
)]
•
(36)
These three isolating integrals are extremely similar to Eqs.(19)-(21): we only have to multiply them with suitable constant factors to obtain the same forms, £,z = - L z ,
(37)
~ ,
(38)
2G
(39)
~_
95
LIU Fu-yao et al. / Chinese Astronomy and Astrophysics 30 (2006) 8~101
Thus, we can transform the 3 isolating integrals Lz, Rx and Ry of the system (27) into the isolating integrals Lz, Rx and Ry of the isotropic linear oscillator (33), and generate the symmetry group SO(3) of rotational transformation in three-dimensional space of the system (33). We now seek a numerical method that will conserve this symmetry group.
3. M A I N T A I N I N G T H E C O N S E R V E D Q U A N T I T I E S OF T H E K E P L E R MOTION WITH THE TRAPEZOIDAL FORMULA
For the first order ordinary differential equation
= f ( t , x),
(40)
one may resort to the linear k-step implicit Newton-Coates quadrature formula, k ~gn+k : 2gn 3- h E
Cikf n+i'
(41)
i=0
where
(_l)k_i ~ok k
Cik -- iV~(k=~) !
I](s - j)ds, j¢i
fn+i = f(tn+i, xn+i),
to obtain numerical solutions, with h the integration step and tn+i = tn 3- ih. Here, a few Newton-Coates integrators are enumerated: h
X n + l -- Xn : - 2 ( f n + f n + l ) ,
k : 1,
]¢ ----2,
k = 3, k=4,
~lgn+3
~n+2 -- ggn ---- h ( f n
--
3- 4 f n + l 3- f n + 2 ) ,
Xn = 3--hs(fn + 3fn+ 1 + 3 f n + 2 3- f n + 3 ) ,
Xn+4 --X~ = -~2h(T f n 3- 3 2 f n + l 3- 12fn+2 3- 3 2 f n + 3 3- 7 fnJc4 ) .
(42)
(43)
(44) (45)
They are called, in order, the second order single-step trapezoidal formula, the fourth order two-step Simpson formula, the fourth order three-step Newton formula and the sixth order four-step Cotes formula. According to Feng Kangls theory [15], it is not difficult to judge that the linear kstep implicit Newton-Coates quadrature formula (41) is symplectic. Feng Kang/s theory is described as follows: For the linear k-step algorithm, k
k
E O~iXn+i = h E ~ifn+i i=O
i=O
(46)
96
LIU Fu-yao et al. / Chinese Astronomy and Astrophysics 30 (2006) 8~101
used to solve linear Hamiltonian systems to be symplectic, it is necessary and sufficient that the two characteristic polynomials of Eq.(46), k
p(~) = ~-~' ai~ ~,
k
a(~) = ~-~.~
i=0
(47)
i=0
should satisfy the relation ¢(~) = - ¢ ( 1 / ~ ) (~(~) = p(~)/a(~)), or that p(~) is antisymmetric while a(~) is symmetric. Therefore we know that Eq.(41) is a linear multi-step symplectic algorithm. What is worth noting is that the symplectic property shown by the linear multi-step formula (41) is tenable only when the method is numerically stable and when the solutions of the difference equation of the linear multi-step method are determined solely by the principal root, free of any effect from the spurious roots [16,17]. Excepting the trapezoidal formula (42) which has a comparatively large and relatively stable region, all other NewtonCoates integration formulae have zero stability, i.e. the stable region is only limited to certain segments symmetrical about the origin on the imaginary axis of the complex plane. For linear systems, so long as the step length is proper, the stability of these linear multistep methods is guaranteed. When they are applied to non-linear systems, rather strong numerical instability will emerge. 3.1 C o n s e r v a t i o n of t h e s y m m e t r y g r o u p of t h e i s o t r o p i c t w o - d i m e n s i o n a l barmonic oscillator The second order single-step trapezoidal formula (42) has a much better stability than other Newton-Coates integration methods and since it is also a sing-step method, it has obvious practical values. But in practice, it has not been given much attention, because it is an implicit method, which generally requires iteration, and so has a lower efficiency than an explicit algorithm. However, when applied to a linear systems an implicit method can be converted into an explicit one. First, for a brief presentation of the problem we shall analyze the harmonic oscillator formula (3) (omitting the subscript i). The differential equation of motion here is (~ = p,
~5= _w2q.
(48)
Its analytic solution is q(t) = qo cos(wt) + ( 1 / w ) p o sin(wt), p(t)
= - w q o sin(wt) + P0 cos(wt).
(49)
Note that Eqs.(49) are only formally different from Eqs. (8) and (9). Applying the trapezoidal formula (42) on Eq.(48), we obtain numerical solutions of the form, qn+l = qn cosc~ + (1/W)pn sin(~, P~+I = -Wqn sin (~ + pn cos c~. Here cosa = ( 4 - w2h2)/(4 +w2h2),
(50)
LIU Fu-yao et al. / Chinese Astronomy and Astrophysics 30 (2006) 87-101
97
sinc~ = 4wh/(4 +w2h2),
~ hf~(h) f~(h) = 4w/(4 - w2h 2) ~-, w + w3h2/4
(51)
It is easily verified that dqn+l A d p n + l = dqn A d p ~ , i.e. Eq.(50) is area-preserving or is a symplectic mapping t. The numerical solution Eq.(50) and the analytic solution Eq.(49) are completely identical in form, but the m a p p e d frequency ft(h) is evidently not equal to the proper frequency w of the system. Therefore, the numerical solution will manifest a linearly increasing along-track error but will not introduce any artificial dissipation. It is not difficult to get H(qn+~,pn+~) - H(q~,pn), showing that the trapezoidal formula strictly conserves the energy of the harmonic oscillator, or produces no truncation error in the energy. This is different from the conservation of energy by the usual symplectic algorithm, which always contains the truncation error in energy which, however, always fluctuates about the original energy, thereby maintaining the energy in a sense, and it m a y be said that the ordinary symplectic algorithm maintains an energy integral in an altered form. Similarly, let us solve numerically the isotropic two-dimensional harmonic oscillator (33) (frequency w = 8 v ~ ) by roans of the trapezoidal formula (42). Then from the result obtained we easily deduce F - 0, and t h a t the 3 isolating integrals Lz, Rx and /~y in Eqs.(37)--(39) are strictly conserved, i.e., the s y m m e t r y group SO(3) they generate is strictly conserved. Then, the isolating integrals in Eqs.(27)--(29) are strictly conserved by the trapezoidal formula (42). 3.2
Conservation of the 5 Kepler orbital elements For the Kepler two-body motion in three-dimensional space, we start with a total of 7 isolating integrals i.e., 1 energy integral, 3 components of the angular m o m e n t u m integral and 3 components of the Runge-Lenz vector (Eqs.(24)-(26)). But since R . L = 0, R and the radial unit v e c t o r / ~ of the periastron point satisfies the relation z
r
--
1 (p.p)r_(r.p)p_#r e#
~ -e
--
=R/(e#)'
--
e#
r - ~ -e
--
--
e#
p= (52)
the Runge-Lenz vector is parallel to the periastron vector and its modulus satisfies the relation IRI = e , (e = IRI/p). (53) Hence there are only 5 independent isolating integrals. (If we use the canonical Delaunay variables, then there will be 5 cyclic coordinates in the new Hamiltonian, producing 5 cyclic integrals and making the existence of 5 independent isolating integrals much more obvious). It is rather obvious t h a t to specify an ellipse centered at the origin in three-dimensional space we need 5 independent parameters, namely, the semi-major radius a, eccentricity e, orbital inclination i, longitude of node f~ and periastron angle w, and hence 5 independent isolating integrals. 4 of the 5 orbital elements other t h a n e can be expressed in terms of the t Note that there exists still the truncation error in the numerical solution which shows up in the mapped frequency ft(h).
98
LIU Fu-yao et al. / Chinese Astronomy and Astrophysics 30 (2006) 87-101
integrals (24)-(26) as follows[26]: a = -#/(2H),
(54)
i = arccos R z ,
(55)
12 = arctan(- Rx / Ry ) ,
(56)
w = arcsin(/3z/sin i) = arcsin[Rz/(e# sin i)].
(57)
Of course, to place ~ and w in their right quadrants, we need the signs of certain trigonometric functions. Minesaki et al. [22] adopted the KS transformation to transform the non-linear Kepler two-body system into the four-dimensional isotropic harmonic oscillator, on which they then applied the trapezoidal formula; they then recovered the original system by inversion, all the while strictly conserving the energy integral, angular momentum integral and Runge-Lenz vector. From this we can conclude that the trapezoidal formula can strictly conserve the 5 Kepler orbital elements a, e, i, ~ and w in the light of Eqs.(53)-(57). Here, we should emphasize that the 5 orbital elements only specify the ellipse in relation to the coordinate origin in three-dimensional space; to specify a particular point on the ellipse, we need a further parameter, the mean anomaly M, say. Now, let us see the effect of the error in mean anomaly M caused by the trapezoidal formula. In Reference [22] a fixed step size is adopted for the physical time t; it is changed to a variable step size for the new coordinate time after the KS transformation. Since a is strictly conserved by the trapezoidal formula, the mean motion n is also strictly conserved. For the physical time covered in one step, the increment of M is A M ----n. At, showing that no truncation error on M is produced by the trapezoidal formula. In contrast, for the usual symplectic algorithm, because there exists a truncation error A n in n, we have ~ = n + A n , and the increment of M is A M = ~ . A t = n . A t + A n . A t , where An. At is the local truncation error in M produced by the usual symplectic algorithm. This shows that the error in M caused by the usual symplectic algorithm induces a longterm variation. This is another aspect where the trapezoidal formula is superior to the usual symplectic algorithm.
4. C O U P L E D
HARMONIC
OSCILLATOR
For the anisotropic linear oscillator, although the original system has the proper frequencies Wl/W2 = n / m , the mapped frequencies (Eq.(51)) obtained through the trapezoidal formula have ~tl(h)/~2(h) ~ n / m . Thus, the trapezoidal formula does not conserve other isolating integrals than the energy integral (or some integral associated with it), that is, the isolating integrals (22) and (23)). Let us take up the coupled harmonic oscillator once again. The coupled harmonic oscillator is H = ~ l ( p ~ +p22 ) + 2lw2r o(.ql2 +q22) - 6qlq2
(58)
with ~ the coupling parameter. Make the coordinate transformation
ql-
q2-
(59)
LIU Fu-yao et al. / Chinese Astronomy and Astrophysics 30 (2006) 87-101
99
and obtain the new Hamiltonian system
1 2 +p~) + ! {022~2 +w~/2) 2 ~ 1~ ,
H = ~(p~
(60)
021 ~- W0v / f l - £/0.)02 ,
(61)
w2 = wo V/I + e/w~ .
(62)
Eq.(60) stands for the anisotropic linear oscillator. When c/w~ << 1, Wl ~ w0[1-e/(2w02)] and ~2 ~ w0[1 + c/(2w0z)]. This indicates t h a t no m a t t e r how small e/w~ is, within an accuracy greater t h a n the accuracy of the computer, the trapezoidal formula no longer conserves the isolating integrals of the coupled harmonic oscillator except the energy integral or any integral associated with it. Thus, the condition of strict conservation of the relevant isolating integrals by the trapezoidal formula is very stringent. Even if the system is a linear system, as soon as it is anisotropic, then the trapezoidal formula ceases to conserve strictly all the isolating integrals except the energy integral or any integral related to it. The Kepler two-body system excepted, all other nonlinear integrable autonomous Hamiltonian systems with n degrees of freedom (such as the one Consisting of a central gravitational field and an additional uniform force field [1] or the problem of two fixed centers) can not, in general, be transformed into an isotropic harmonic oscillator. If there exist more t h a n n independent isolating integrals, the trapezoidal formula generally can strictly conserve at most n such integrals.
5. C O N C L U D I N G
REMARKS
The linear autonomous Hamiltonian system with n degrees of freedom certainly has n mutually involutive independent isolating integrals, and these independent isolating integrals are strictly conserved by the trapezoidal formula. If there exists commensurability among some or all of the characteristic frequencies, then the number of independent isolating integrals can exceed n, but can not exceed 2n - 1. When the system is also isotropic, it possesses 2 n - 1 independent isolating integrals, and all of t h e m can be strictly conserved by the trapezoidal formula. If all the characteristic frequencies of the system are commensurate and the system is isotropic then even though there exist 2n - 1 independent isolating integrals, none of these are conserved by the trapezoidal formula except the energy integral (or some other integral related to it). The two-dimensional Kepler two-body system pertains to an isotropic nonlinear system and has the energy integral, the angular m o m e n t u m integral and the Runge-Lenz vector. In the case of plane motion, among these three quantities there are a total of 4 nonzero isolating integrals but only 3 of them are independent ones. These 3 independent isolating integrals make up the s y m m e t r y group SO(3) of the system, and via the Levi-Civita transformation it cam be transformed into the s y m m e t r y group of the two-dimensional isotropic harmonic oscillator, which is strictly conserved by the trapezoidal formula. In the case of motion in three-dimensional space, the three n a m e d quantities have in total 7 isolating integrals of which only 5 are independent, as the Runge-Lenz vector is parallel to the the periastron
100
LIU Fu-yao et al. / Chinese Astronomy and Astrophysics 30 (2006) 8"7101
point vector. Through the KS transformation the system can be transformed into the four-dimensional isotropic harmonic oscillator (in topological language, the two- and threedimensional Kepler two-body systems are respectively homeomorphic to the two- and fourdimensional isotropic harmonic oscillators). The trapezoidal formula can likewise strictly conserve the three named quantities (the energy and angular momentum integrals and the Runge-Lenz vector), and this fact is reflected, in particular, in the 5 Kepler orbital elements a, e, i, gt and w. As for other nonlinear integrable autonomous Hamiltonian systems, generally speaking, they can not be transformed into isotropic harmonic oscillators, and the trapezoidal formula can not in general conserve isolating integrals with forms similar to the form of the Runge-Lenz vector. It needs to be emphatically pointed out that although it is seen analytically that the trapezoidal formula strictly conserves the energy integral, angular momentum integral and Runge-Lenz vector as well as the 5 Kepler orbital elements a, e, i, ~t and w, numerical results obtained by the computer always have errors that increase linearly with the integration time. This is because freedom from the truncation error makes the rounding error the principal source of error. In passing, it may be mentioned that the equations of stellar orbital motion and of the relativistic perihelion precession can be simply and conveniently derived by taking advantage of the Runge-Lenz vector [27]. ACKNOWLEDGEMENT The second author of this article discussed the concepts of isolating and non-isolating integrals with Prof. H U A N G Tian-yi of Department of Astronomy of Nanjing University and derived great benefits from the discussion. We wish to express our indebtedness to Prof. Huang. References 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18
Arnold V.I., Mathematical Methods of Classical Mechanics (translated by Vogtmann K. and Weinstein A.). New York: Springer-Verlag, 1978, 100 LichtenbergS.J. Lieberman M.A., Regular and Chaotic Dynamics. Springer-Verlag, 1983, 1 GoldsteinH., Classical Mechanics. Addison-Wesley Publishing Company, 1950, 55-300 Wu X., Huang T.Y, Wan X.S., ChASzA, 2003, 27(1). 114 Wu X., Huang T.Y., Wan X.S., AcASn, 2002, 43(4), 391 Liu F.Y., Wu X., Lu B.K., ChA&A, 2005, 29(1), 92 Liu F.Y., Wu X., Lu B.K., AcASn, 2004, 45(4), 402 Wu X., Huang T.Y., Wan X.S., ChA&A, 2005, 29(1), 81 Wu X., Huang T.Y., Wan X.S., AcASn, 2004, 45(3), 310 WisdomJ., Holman M., AJ, 1991, 102, 1528 YoshidaH., CeMDA, 1993, 56, 27 AscherU.M., Numerical Algorithms, 1997, 14, 1 YoshidaH., Phys. Lett. A, 2001, 282, 276 YoshidaH., CeMDA, 2002, 83, 355 Feng K., Proceedings of the Annual Meeting on Computational Mathematics in Tianjin. Tianjin, Nankai University Press, 1990, 1 Huang T.Y., Wang C.B., Zhao Z.Y., ChA&:A, 1998, 22, 101 Huang T.Y., Wang C.B., Zhao Z.Y., AcASn, 1997, 38, 278 Chiou J.C., Wu S.D., J. Chem. Phys., 1997, 107(17), 6894
LIU Fu-yao et al. / Chinese Astronomy and Astrophysics 30 (2006) 87-101
19 20 21 22 23 24 25 26 27
101
Zhu W., Zhao X., Tang Y., J. Chem. Phys., 1996, 104(6), 2275 Kalogiratou Z., Simos T.E., J. Comput. ~ App. Math., 2003, 158, 75 Minesaki Y., Nakamura Y., Phys. Lett. A., 2002,306, 127 Minesaki Y., Nakamura Y., Phys. Lett. A., 2004, 324, 282 Hdnon M., Heiles C., AJ, 1964, 69, 73 http://wwwgoogle.com. Annu. Rev. Astron. Astrophys., 1982, 20, 399: Casati G., Chirikov B.V., Ford J., Phys. Lett. A., 1980, 77, 91 Brouwer D., Clemence G.M., Methods of Celestial Mechanics, New York: Academic Press, 1961, 15 Weinberg S., Gravitation and Cosmology, New York: John Wiley, 1972, 165