A pseudo third order symplectic integrator

CHINESE ASTRONOMY AND ASTROPHYSICS ELSEVIER

Chinese Astronomy and Astrophysics 29 (2005) 92-103

A Pseudo

Third Order Symplectic Integrator t *

L I U F u - y a o l'a

WU Xin 2

L U B e n - k u i 1,3

1purple Mountain Observatory, Chinese Academy of Sciences, Nanjing 210008 2School of Sciences, Nanchan 9 University, Nanchang 330047 3National Astronomical Observatories, Chinese Academy of Sciences, Beijing 100012

A b s t r a c t The symplectic integrator has been regarded as one of the optimal tools for research on qualitative secular evolution of Hamiltonian systems in solar system dynamics. An integrable and separate Hamiltonian system H = Ho + Y]N=1 eiHi (e/ << 1) forms a pseudo third order symplectic integrator, whose accuracy is approximately equal to that of the first order corrector of the Wisdom-Holman second order symplectic integrator or that of the ForestRuth fourth order symplectic integrator. In addition, the symplectic algorithm with force gradients is also suited to the treatment of the Hamfltonian system H = Ho(q, p) + ell1 (q), with accuracy better than that of the original symplectic integrator but not superior to that of the corresponding pseudo higher order symplectic integrator. K e y w o r d s : celestial mechanics--methods: numerical

1. I N T R O D U C T I O N The symplectic algorithm is the integration method of finding numerical solutions of Hamiltonian canonical equations while the symplectic structure dq A dp of the Hamiltonian system is maintained. The Newtonian n-body problem in solar system dynamics belongs to the Hamiltonian problem. Therefore, the symplectic algorithm has been widely applied and developed in the qualitative research of celestial mechanics. t Supported by National Natural Science Youth Foundation Received 2004-01-09; revised version 2004-04-07 * A translation of Acta Astron. Sin. Vol. 45, No. 4, pp. 402-412, 2004 0275-1062/05/S-seefront matter © 2005 Elsevier B. V. All rights reserved. DOI: 10.1016/j.chinastron.2005.01.009

LIU Fu-yao et al. / Chinese Astronomy and Astrophysics 29(2005) 92-103

93

The symplectic method is of two types, implicit and explicit, and the latter requires the Hamiltonian function to be of the integrable form. For the Hamiltonian system of the form H = Ho + ell1, the truncation error Herr is a power series in the step length ~- and the small parameter e with multiple Poisson brackets of H0 and/-/1 as coefficients[1]. In Ref. [2] four ways of improving the accuracy of the symplectic integrator are outlined: (1) Raising the power of ~- in Herr, i.e. constructing higher order series, such as is done in Yoshida [3] and Forest et al. [4]. Symplectic integrators above the 4 th order are seldom adopted because of the decrease in computational efficiency when there are too many substeps. This is the basic reason why the accuracy of the symplectic integrator is not high enough. (2) The way the Hamiltonian function is decomposed has a direct effect on the accuracy of the symplectic integrator. The smaller e is, the better the numerical accuracy. In the decomposition H = T + V, there is no difference between the kinetic energy T and the potential energy V in order of magnitude, and Herr of the m-th order symplectic integrator is O(~-rn). Wisdom et al.[5] decomposed the Hamiltonian function of the n-body problem of the solar system into a primary, integrable Kepler part and a secondary, interplanetary interactive part by taking advantage of the Jacobi coordinates, making Herr be O(c~'m), so that for the same step length its error is e order of magnitude smaller than that of the original formulation, with notable improvement in the computational efficiency and the accuracy. Duncan et al. [6] redecomposed the Hamiltonian function of the n-body problem of the solar system into three integrable parts, which has also the merits of similar improvement on the accuracy of the symplecctic integrator. Wu et al. [7] researched the symplectic algorithm of the Hamiltonian system H = H0 + ~N= 1 eiHi(ei << 1), decomposed it into many integrable parts, and found an almost ideal integrator. (3) The integration can be similarly optimized without decreasing the step length or increasing the order of integration. For example, the terms including the first and second powers of c in Herr may be eliminated,-- this is equivalent to raising the power of ~ in Herr, and is the very symplectic correction of Wisdom et al. is]. It does not play a role in the whole process of the numerical integration and only changes the initial condition and the value of the output point, it both maintains the conservation of energy and improves the computational efficiency. Wu Xin et al.[9] extended the symplectic correction to the case where the Hamiltonian function is decomposed into a number of integable parts. When only terms of the first power of e are eliminated, we obtain a pseudo high-order symplectic integrator [1°], composed of a limitd number of symmetrically combined Lie operators. Therefore, when e is small, the pseudo high-order symplectic integrator is greatly superior in efficiency and accuracy to the standard integrator of the same order. (4) The step length ~- has an effect on the accuracy of integration: the smaller ~- is, the more accuracy the integration. However, if 7 is too small, the integration time and hence the cumulative error would be excessive. Thus, for a given integrator, a suitable step length needs to be selected. According to the Wisdom-Holman leap-frog mapping a suitable step length for the outer solar system is 350 days [51. On the other hand, different celestial bodies move with different periods, and different step lengths should be adopted for the different bodies while maintaining the symplectic structure (see Ref.[6]). Also, at a close encounter of a minor planet with a major planet or when the minor planet is at the the pericentre of a very eccentric orbit, the gravitation on the minor planet exerted by the central body is very great, the step of integration needs to be changed promptly, otherwise the numerical solution will not be realistic. This is the technique of normalization [111 that maintains the

94

LIU Fu-yao et al. / Chinese Astronomy and Astrophysics 29 (2005) 92-103

symplectic structure [12-15]. The outstanding advantage of this sort of symplectic algorithm is maintenance of a good numerical stability in case of highly eccentric orbits. Generally speaking, the above-mentioned symplectic integrators can operate with the time reversed. However, for irreversible evolutionary systems such as the quantum mechanics multiple body problem or attractors, we have available only the first and second order symplectic integrators and the pseudo high order symplectic integrators; the other integrators can not be applied. This is because negative coefficients appear in the combinatorial operators in the higher order methods of Yoshida [3] and Forest et al. [4] and in the symplectic correction, or negative step lengths are adopted. For the case of decomposition of H = T + V, Ruth [16] constructed a positive coefficient third order symplectic integrator, and later on the introduction of the gradient of force through the work of Suzuki [17] and Chin [ls'19] et al. made possible a positive step length high order symplectic algorithm. It is pointed out in Refs. [18] and [19] that the accuracy of the positive step length high order symplectic algorithm containing the gradient of force in the case of equal step length is from 10 to 8(} times better than that of the Forest-Ruth method of the same order [4]. When the integrable decompositions of the Hamiltonian differ by orders of magnitude, the pseudo high order symplectic integrator is the best choice because of its accuracy and efficiency and the simplicity of its program. Based on this consideration, a pseudo third order symplectic integrator is constructed using asymmetrical operators in this article. To evaluate the merit of the integrator, we have numerically compared it with the second order symplectic algorithm, the first symplectic corrector and the Forest-Ruth fourth order symplectic method [4]. Naturally, the pseudo high order symplectic integrator should also be compared with the corresponding positive step length symplectic algorithm of the same order and with gradient of force. Thus we explore in this article the feasibility of dealing with the case when the decomposed parts of the Hamiltonian function mutually differ greatly by means of the symplectic integrator with the gradient of force.

2. A P S E U D O T H I R D O R D E R S Y M P L E C T I C I N T E G R A T O R Consider the Hamiltonian system H = Ho + ¢HI and assume that both Ho and HI are integrable and ~ = IHI/Hol ~ i. Let A and B be the Lie derivatives of Ho and //i, A = LHo = {H0}i B = LH1 = {H1} and W = ~-L~ = T{H, f}, {H, f} being the Poisson

brackets of the functions H and f. The forward numerical solution for initial conditions OH

(qo, P0) of the canonical equations q = -~p and/b -

OH

Oq is

where [a] k exp(W) = exp[~-(A + ~B)] ,.~ H exp(bi~-~B) . exp(ai~-A).

(1)

The symplectic integrator of any order can be established by making a choice of the coefficients a~ and b~ When i. = 1, the combinatorial operator is

LIU Fu-yao et al. / Chinese Astronomy and Astrophysics 29(2005) 92-103

exp(TeblB)- exp(~-alA) = exp(W1).

95

(2)

According to the Baker-Campbell-Hausdroff (BCH) formula [2°1, one may have Wl

:

T( al n ~- £D1B ) -~- ~ 72 eal b1[B, A]

+ 1 T 3 ( e2al b~[B, B, A] +

ea2bl [A, A, B]) + l'r4e2a~b2[B, A, A, B] 1 .rh(ea4bl[B, B, B, B, A] + e4alb4[B, B, B, B, A]) + 720 1 .rh(e3a~b~[A' B, B, B, A] + c2a~b2[B, A, A, A, B]) +

(3)

36O

1 Th/ 3 2 - 3 r ~ - ~ 120 Lc alOl[~'B'A'A'B] + ¢2a31b~[A,A,B,B,A] ) + " " .

It is stipulated that the commutators [B, A] -- B A - AB, [B, B, A] = [B, [B, A]] and so forth. If al = bl = 1, then the first term of the above equation just corresponds to the original Hamiltonian while the other terms are the error part of the solution. Evidently, the combinatorial operator (2) is the first order symplectic integrator. For the sake of simple and convenient writing, we find the "product" of the exponential operators by making use of the (BCH) formula and retain only terms up to 0(7-3). Writing exp(Ta2A).exp(W1) = exp(W2), we have

W2 = "5(al + a2)d + ~'~blB + l~-2ebl(al - a2)[B, A] +

--~T3e2(al + a2)bl [B, B, A] +

(4)

~--~-3ebl (a~ - 4ala2 + a~)[A, A, B] + O(cT4). Take al + a2 = 1, bl = 1 and bl(ai - a2) = 0, or al = a2 = 1/2 and bl = 1, and we obtain the symmetric second order symplectic combinatorial operator often used in literature T

T

$2 = e x p ( ~ A ) , exp(TcB) • exp(~A).

(5)

The error part of its solution is 1 3e[A, A, B] + O(CT4), ERRs2 = l~-3e2[S, B, A] - -~7

(6)

When i = 2, let exp(Teb2B), exp(W2) = exp(W3), then one may have W3 ~- ( a l + a2)'rA + (bl +

b2)'reB +

2[(al + a2)(bl + b2) - 2a2bl]T2e[B, A] + 1 [ ( a l + a2)(bl + b2) 2 - 6a25152173e2[B,B, A] + 1 ~ [ ( a l + a2)2(bl + b2) - 6ala2bl]T3c[A, A, B] + O(e74).

(7)

96

LIU Fu-yao et al. / Chinese Astronomy and Astrophysics 29 (2005) 92-103

The expression of W3 has appeared in Ref. [18] . We now consider in detail three cases of the positive coefficients al, a2, bl and b2. 1) W h e n a x + a 2 = 1, b l + b ~ = 1 and2a2bl = 1, i.e. el = 1 - a 2 , 51 = 1/(2a2) and b2 = 1 - 1/(2a2), infinite second order symplectic integrators can be obtained from exp(W3). 2) When al + a2 = 1, bl + b2 = 1, 2a2b1 = 1 and 6a2blb2 = 1, or al -- 1/4, a2 = 3/4, bl = 2/3 and b2 -- 1/3, exp(W3) can be written as S~ = exP(3~'eB ) • e x p ( ~ - A ) , exp(~TcB) . 2 exP(4TA) .

(8)

The error of its solution is ERRz~ = 1T3c[A, A, B] + O ( e T 4 ) . (9) 48 This shows that S~ still belongs to the second order algorithm. 3) When al + a2 = 1, bl + b2 = 1, 2a2bl = 1 and 6ala261 = 1, or al -- 1/3, a2 = 2/3, bl = 3/4 and b2 = 1/4, exp(W3) has the form

PS3

=

exp(17-eB) • e x p ( h2 T d ) . exp(3 eB) exp(hTA 1 ).

(10)

The error of its solution is ERRps3 = ~8T3e2[B, B, A] + 0(c74).

(11)

In other words, the error of PS3 is max{O(~-3e2), O(e~-4)}. If the error is O(c~-4), then PS3 becomes a third order integrator; if the error i s O(T3e2), it is by e orders of magnitude smaller than Eq.(9), also equivalent to a raising the of order of the integrator. Essentially PS3 is part of the second order algorithm. However, since c is smaller and its accuracy is obviously better than that of $2, PS3 has the accuracy of the third order integrator. Thus, PS3 is called a "pseudo third order symplectic integrator". Of course, various higher order symplectic integrators can be similarly formed by using the asymmetrical combinatorial operator (1), but the derivation is rather complicated. Thus, in most of the references the symmetrical operators are adopted and the terms of even powers of 7 are not included in the error ERR. The examples are: Yoshida [3]'s high order sympleetic method, Forest-Ruth [4]'s high order method and Chambers-MurisonD°]'s symplectic integrators with orders higher than the third. Here, we cite the Forest-Ruth fourth order symplectic method: $4 = e x p ( ~ c B ) . e x p ( ~ A ) , e x p ( T c k T e B ) . e x p ( - k ~ - A ) -

where k = 21/3 and c = 2 - k. In addition, a Chambers-Murison pseudo fourth order symplectic integrator is given as follows: T£

T

2T£

T

T~

PS4 = exp(-~-B) - exp(~A) • exp(--3--B ) • e x p ( ~ d ) • exp(--~-B) T3

~T 5

= exp[T(A + c B ) ÷ --~e2[B,B,A] + -~6[A,A,A,A,B] ÷...].

(13)

LIU Fu-yao et al. / Chinese Astronomy and Astrophysics 29(2005) 92-103

97

In the next section we shall evaluate and estimate the numerical effects of the pseudo third order symplectie integrator PS3, the second order symplectic algorithm $2 and its first symplectic correctors C~q2, S~, $4 and PS4 by means of numerical experiment.

3. N U M E R I C A L

EXPERIMENT

In this section the method presented in Ref. [9] is adopted, the Sun-Jupiter-Saturn three body problem is taken as the model, and the various above-mentioned algorithms are compared from numerical simulation regarding the energy error and tracking error, the latter being one of the most important measures of the accuracy of the integration. For convenience, the tracking error is replaced by the mean longitude error. In terms of the heliocentric coordinates Qi and the barycentric momentum P i (i : 0 stands for the Sun), Duncan et al.[6] decomposed the Hamiltonian function in the barycentric coordinates for the solar and planetary n-body system as: H = Hkep + Hsun -~ H i n t ,

(14)

where n--1

V'(IIP II:2

Hkep = z.._~ i=l

2mi

amain0 ItQill )'

(15)

n--1

/-/sun

1

-- 2m011 ~P~l[ 2, i=l

n - 2 n-1

Hint i=1

j=i-t-1

(16)

Gmimj

IJ -- ll"

(17)

Here mi is the mass of the i-th body. Hkep stands for the independent two-body (Sun and planet) motion, of which the solution is integrable and can be expressed by the Gauss function. //sun is the part of the motion of the Sun with respect to the barycentre, which is integrable because it includes only the generalized momentum. H i n t is the interactive potential between the planets, which is a function of the coordinates only and is therefore also integrable. As for the Sun-Jupiter-Saturn three body system, since the four interior planets and the Sun represent an approximately symmetrical distribution, the masses of the four interior planets are simply added to the Sun. Here, eI = IHsun/Hkep] and e2 = [Hint/Ykepl are each about 10 -4. Denote A = {Hkep}, B1 = {//sun} and B2 = {Hint}- The above exponential operator exp(~-eB)[see Eqs.(5), (8), (10), (12)and (13)] is changed into exp(~'e2B2).exp(~-el B1). The DE405 ephemeris heliocentric coordinates are taken as the initial condition and the computed results using the high-accuracy Cowell 12th order method are used for the reference orbit. The step length in the Cowell method is taken as 35 days, while for all the sympleetie algorithms, the step length is fixed at 350 days. The integration is taken to 10 s years and the data is stored about every 20,000 years. The mean longitude error AA of the planet is taken to be the difference between the calculated results of the Cowell method and of the

98

LIU Fu-yao et al. / Chinese Astronomy and Astrophysics 29 (2005) 92-103

symplectic algorithm. It should be pointed out that the position coordinates of the Sun with respect to the barycentre of the solar system can not be directly found from the numerical integration, but can be obtained by making use of the conservation of momentum. Fig.1 shows the relative error in the energy A E / E o , as a function of time, for the six symplectic algorithms. The fact that, in none of them is there a secular variation fully demonstrates the important characteristic of the symplectic algorithm. The energy is best maintained by the first symplectic corrector, CS2, of the second order symplectic integrator $2, with an error on the order of 10 -7. The Chambers-Murison pseudo fourth order symplectic integrator, PS4, has the same order of error. For the other algorithms, the errors are all of the order of 10 -6, with the pseudo third order symplectic integrator PS3 slightly better than the Forest-l~uth fourth order symplectic method $4. As is pointed out above, S~ is altogether a second order algorithm, and its energy error is slightly smaller than that of $2. The reason is that the absolute value of the coefficient of the commutator [A, A, B] in Eq.(9) is one-half that in Eq.(6). PS3 is nominally a second order algorithm but it shows the accuracy of a third order method because the error in Eq.(11) contains the term of the second power of the small disturbing parameter e = max{el, e2}. It is indicated that PS3 is one order of magnitude inferior to CS2, but slightly better than $4 in regad to the energy error. However, the along-track error is not completely so. The time evolution of the mean longitude error AA of Jupiter in the integration by the six symplectic integrators is shown in log-log plots in Fig.2. AA is restricted to the interval [0, 2~]. Before AA reaches the first cusp~ it increases almost linearly with time in all the six figures with no exception. This is also one of the important characteristics of the symplectic algorithm. The time when the mean longitude error AA first reaches 2~r is denoted by T, which can serve as measure of how good the integrator is (see Refs.[2,7,9]): the larger T is, the more slowly the mean longitude error AA increases, and the more accurate the integrator is. Table 1 lists the values of T for the six integrators. It shows that the numerical stability of PS3 is one order of magnitude better than that of $2 while being of the same order as those of CS2 and $4, with little difference between them. As far as computational efficiency (the CPU times taken by the various integrators are listed in Table 2) is concerned, PS3 is similar to the second order symplectic algorithms $2 and S~. Table 1

T h e t i m e T (unit: 106 years) taken by various s y m p l e c t i c i n t e g r a t o r s w h e n the m e a n longitude e r r o r AA of the J u p i t e r reaches 27r algorithm T

Table 2

$2 3.3

CS2 69

S~ 5.7

PS3 52

$4 74

PS4 > 100

T h e C P U times (hour-" minute: second) consumed by various s y m p l e c t i c algorithms algorithm CPU time

$2 0:46:15

CS2 0:54:12

S~ 0:48:09

PS3 0:48:07

$4 1:48:33

PS4 0:51:11

99

LIU Fu-yao et al. / Chinese Astronomy and Astrophysics 29(P005) 92-103

~f~

~. 2 v

"~ -2 -4 4

6

8

10

0

2

t(I 0 7 years)

4

6

8

10

8

10

8

10

t(10 7 years) 1.0

,_, 0.5 0 rrl "~ -1

"g 0.0

-0.5 -2 4

6

8

10

0

2

t(l 0 7 years)

4

6

t(lOTyears) 1.O 0.5 '7

0.0 -0.5 -1.0

rT1

-1.5 -2.0 --v 0

4

6

t(10 7 years)

8

10

-2.5

0

2

4

6

t(10 7 vears/

Fig. 1 T h e relative e n e r g y error AE/Eo as a function of t i m e in t h e i n t e g r a t i o n of t h e S u n J u p i t e r - S a t u r n t h r e e b o d y s y s t e m by various s y m p l e c t i c integrators. It c a n be seen t h a t t h e p s e u d o t h i r d order s y m p l e c t i c i n t e g r a t o r PS3 is slightly s u p e r i o r to t h e F o r e s t - R u t h f o u r t h order s y m p l e c t i c i n t e g r a t o r $4, b u t g r e a t l y inferior to t h e first order corrector CS2 of t h e second order s y m p l e c t i c i n t e g r a t o r S2.

LIU Fu-yao et al. / Chinese Astronomy and Astrophysics 29 (2005) 92-103

100

.<] o

1.0 0.5 0.0 -0.5 -1.0 -1.5 -2.0 -2.5 -3.0

/4.5

1.0 0.5 0.0 -0.5 -1,0 -1.5 -2.0 -2.5 -3.0

j "

5,0 5.5

6.0 6,5 7.0 7.5

8.0

4.5

5,0 5.5 6.0 6,5 7.0 7.5

Ioglo t

loglo t 1.0 0.5 0.0 -0.5

1,0 0.5 0.0 -0.5 <] o -1,0 -1,5

% -1.o

-2.0

-1.5 -2.0

-2.5 -3.0

-2.5 -3.0

4,5

5,0 5.5

6.0 6.5 7.0 7.5

8.0

4.5

5.0 5.5

loglo t

6.0 6.5 7.0 7.5

8.0

loglo t LO 0.5 0.0 -0.5

1.0 0.5 0.0 -0.5 % -1,0 -1.5 -2.0 -2.5 -3.0

8.0

% -1.o -I.5 -

4.5

5.0 5.5

6.0 6.5 7.0 7.5 loglo t

8.0

-2.C

-2.5 -3,C -3.5

4.5

5.0 5.5

6.0 6.5 7.0 7.5 loglo t

8,0

Fig. 2 The time curve of the mean longitude error AA (denoted by A) of Jupiter produced when the Sun-Jupiter-Saturn three body system is integrated with the various symplectic integrators. It can be seen t h a t the accuracy of the pseudo third order symplectic integrator PS3 is about equal to t h a t of the Forest-Ruth fourth order symplectic integrator ,5'4 or the first order corrector CS2 of the second order symplectic integrator $2.

L I U Fu-yao et hi.

4. C O M P A R I S O N

/

C h i n e s e A s t r o n o m y and Astrophysics 29(2005) 9 2 - 1 0 3

OF PSa W I T H T H E S T A N D A R D S Y M P L E C T I C M E T H O D Sa

I01

THIRD ORDER

In the preceding sections the accuracy of PS3 has mainly been compared with the accuracies of the second order symplectic algorithm $2 and its first order symplectic corrector CS2, the fourth order method $4 and the pseudo fourth order symplectic integrator PS4. Of course, the accuracy of PS3 should naturally be compared with that of the standard third order method $3, which is not in the form of a symmetrical combinatorial operator. Sa can be conveniently constructed using the operator combinatorial expression of PS3. Note, in Eq.(11) for the error in the solution of PS3, [B, A, B] has as coefficient -1/48. Therefore, in Eq.(10) we need only substitute/) for B in the operator exP(4~-e.B ) and have

= eB + 17-2e2[B,A,B].

(18)

This is precisely the positive coefficient third order symptectic integrator established by Ruth[161 $3 = exp(17-B) • exp(g~-A) 2 1 • exp(~7-eB) • exp(~7-A).

(19)

However, Ruth adopted the form of H = T(p) + V(q) for the decomposition of the Hamiltonian. The Lie operators corresponding to T(p) and V(q) are written as T and V, respectively. Substituting T and V for A and eB in Eq,(19), we have i

~

2

3

1

TVS3 = exp(~7-V) • exp(gT-T ) • exp(47-V) • exp(~7-T) ,

(20)

where

= V + 1--7-2[V,T,V].

(21) 12 Here lies the difficulty of the actual computation of the two-fold operator [I7,T, V], and this is essentially the reason why this kind of symplectic integrator has received scant attention. The gradient of force introduced by Chin [ls,19] makes it possible to realize the positive step length higher-order symplectic algorithm. The method is as follows:

0 _ ~ V~llfll2 0 IV, T, V] = 2 ~ L OL Oq~ 0;~ I----"M i

i,j

(22)

~p i

where fj(q) = -OV(q)/Oqj. The difference scheme of the third order symplectic combinatorial operator TVS3 is 7-

q* = qn-1 + "~Pn-1,

p* = Pn-I + ~ Y ( q * ) , 2T

qn = q* + ' T P 7-

.

(23)

, 2

Pn = P* + ~($(qn) + ~2VllY(q~)Jl2).

102

LIU Fu-yao et al. / Chinese Astronomy and Astrophysics 29 (2005) 92-103

Besides, Chin [18,19] has also constructed several fourth order symplectic integrators with the coefficients or step lengths all positive. Unlike the negative coefficient higher order methods of Yoshida [31 and Forest [4], they can deal with the irreversible evolutionary systems. When the Hamiltonian is decomposed in the Wisdom-Holman[ 5] form, H = Ho(q, p) + eHl(q), we believe that the third order symplectic combinatorial operator $3 can also be realized effectively. We still assume fj (q) = -OeHl(q)/Oqj. The Kepler part Ho(q, p) in the time ~- undergone by the initial state (q0, P0) has the analytic solution, q = Kepler(q0, P0, ~-) and p = Kepler(qo,Po , "r). In order to find the solution o f / ) , the operator A in Eq.(18) is decomposed into T + V, where the Hamiltonian functions corresponding to T and V are the kinetic energy part and potential energy part of H0. Then

c2[B,A,B]

e2[B, T + V, B] = e2 [B, T, B] + e2 [B, 17, B] = ofj

o

_

2 E.. fJ Oqi Opi ~,3

Zv llfll i

0

Opt'

(24)

where both the operators V and B in [B, V, B] are used only for the solutions of the momentum p. Thereby V B = BV, or [B, V, B] is the null operator. The symplectic difference scheme of Sa is T

q* = Kepler(qn_l,pn_l, ~),

p* = Kepler(qn_l,pn_l, ~), p ** = p* + "J--~f (q *) , q~ = Kepler(q*,p**, ~2r- ) ,

(25)

27 p*** = Kepler(q*,p**,-3-),

p~ = p*** + "~(f(q~) +

vllf(q~)ll2).

Evidently, the accuracy of $3 is much better than that of TVS3 for the same step length. We note that the error of PSa is max{O(e2~-a), O(eZ4)}, hence a proper choice of e and can make the error of PS3 to be O(er4). Thus, Sa is not superior in comparison with PS3 in accuracy and it is rather complicated to find the gradient of force for $3 when a multiple body problem is discussed. By the way, this approach can also be applied to the fourth order symplectic integrator with the gradient of force in Refs.[18,19].

5. C O N C L U S I O N S For the n-body problem of the solar system, if we adopt either the Wisdom-Holman or the Duncan et al. decomposition of the Hamiltonian function, then the pseudo third order symplectic integrator has the advantages of involving fewer substeps, simpler programming and higher computational efficiency, and especially, it reaches the accuracy of the standard

LIU Fu-yao et al. / Chinese Astronomy and Astrophysics 29(2005) 92-103

103

third order symplectic integrator, that is, about the same accuracy as the first order corrector of the Wisdom-Holman second order symplectic integrator or the Forest-Ruth fourth order symploetic algorithm. Thus, PS3 has practical values and is hereby recommended. References 1 2 3 4 5 6 7 8 9 10 11

12 13 14 15 16 17 18 19 20

Yoshida H., CeMDA, 1993, 56, 27 Wu X., Symplectic Algorithm and the Chaotic Phenomenon in the Geneeral Relativity, a doctoral dissertation of Nanjing Univeristy, 2003 Yoshida H., Phys. Lett.A, 1990, 150, 262 Forest E., Ruth R.D., Physica D, 1990, 43, 105 Wisdom J., Holman M., A J, 1991, 102, 1528 Duncan M.J., Levision H.F., Lee M.H., AJ, 1998, 116, 2067 Wu X., Huang T.Y., Zhang H., et al., Ap&:SS, 2003, 283(1), 53 Wisdom J., Holman M., Touma J., Fields Inst. Commun., 1996, 10, 217 Wu X., Huang T.Y., Wan X.S., AcASn, 2002, 43(4), 391 Chambers J.E., Murision M.A., AJ, 2000, 119, 425 Baumgarte J., Stiefel E., Examples of Transformations Improving Numerical Accuracy of the Integration of Differential Equations, in Proceedings of the Conference on the Numerical Solution of Ordinary Differential Equations, Dale ed., Springer-Verlag, 1974 Mikkola S., CeMDA, 1997, 68, 249 Mikkola S.M., Aarseth S., CeMDA, 2002, 84, 343 Preto M., Tremaine S., A J, 1999, 118, 1532 Rauch K.P., Holman M., AJ, 1999, 117, 1087 Ruth R.D., IEEE Trans Nucl Sci, 1983, NS-30, 2669 Suzuki M., Phys. Lett.A, 1995, 201,425 Chin S.A., Phys. Lett.A, 1997226344 Chin S.A., Chen C.R., arXiv: astro-ph/0304223, 2003 Vaxadarajan V.S., Lie Groups, Lie Algebra and Their Representation, Englewood Cliffs, 1974