LITERAL

SOLUTION

FOR

DIETER

HILL'S

LUNAR

PROBLEM

S. S C H M I D T

University of Cincinnati, Department of Mathematical Sciences, Cincinnati, Ohio 45221, U.S.A.

(Accepted 20 October, 1979)

Abstract. A computer p r o g r a m for the manipulation of power series in several variables is used to find the first order solution to Hill's lunar problem. The ratio m of the mean motion of the Sun to that o f the M o o n is kept as a formal parameter. The series in m are k n o w n to converge very poorly. It is shown how efficient algorithms among them the Lie transformation allow us to c o m p u t e the series in m as far as they are needed. W h e n the series are evaluated at Brown's numerical value for m the results achieve or exceed his accuracy.

1. Introduction

The main problem of lunar theory consists of finding the solution to a three body problem under the assumption that the distance between two bodies (Earth and Moon) remains much smaller than their respective distances to the third body (Sun). Brown (1899-1908) gave an approximate solution to this problem in the form of a power series in several parameters. His work is based on two fundamental papers of Hill (1878, 1886) who had considered a slightly simplified problem by assuming that the Sun and the Earth-Moon system rotate around each other uniformly and that the Sun is infinitely far away. One of the parameters in the lunar theory is the ratio of the mean motion n' of the Sun to the mean motion n of the Moon in a rotating coordinate system. It is denoted by m

=

n'/(n

-

n').

Series expansions in this parameter converge very poorly. Since m is known to a fair degree of accuracy and to avoid the problem of convergence, Brown and Hill from the beginning used the numerical value m = 0.08084 89338 08312

(1)

whenever their computations had to be accurate for practical applications. Hill (1894) gave the series expansion in m for the motion of the lunar perigee through m 9 and through eleventh order in a related parameter. He had to conclude that the series expansion was not far enough to give results which were comparable to the method of special value. Bourne (1972) repeated Hill's computation and found some numerical errors but he did not push the series any further. We will show how to find the series in m for a complete first order solution of Hill's lunar problem. Since the computations will be done by machine, methods will be used Celestial Mechanics 19 (1979) 279-289. 0008-8714/79/0193-0279501.65. Copyright 9 1979 by D. Reidel Publishing Co., Dordrecht, Holland, and Boston U.S.A.

280

D I E T E R S. SCH_IVLIDT

which are different from those employed by Hill. In particular the method of Lie transformation will be used to solve Hill's differential equation. Efficient algorithms enable us to compute the series as far as necessary so that we obtain the same result as Brown and Hill when we replace m by its numerical value. For example, for the motion of the perigee we find terms through m a~ and the solution for the variational equations is given through m 29. Cowell (1896) gave the solution for the inclination but in literal form only through sixth order in m. It is no effort to extend his series to order sixteen by machine so that we get the accuracy which Brown achieved in his lunar theory. All computations were performed on the Amdahl 470-V6 at the University of Cincinnati. The machine is running under the IBM operating system OS/VS2 release 1.7. The program is written in the computer language PL/I. It is called as a subroutine from a package of PL/I programs developed for the manipulation of real or complex power series in several variables. The package is named P O L Y P A K and was derived from lectures given by A. Deprit at the University of Cincinnati on an automated series processor. As a reference see Rom (1969, 1971). The latest version of Polypak uses balanced binary trees for the insertion of terms into a series see Knuth (1968). Total computing time is about 12 min for finding the intermediate orbit through m 18, the solution for the inclination through m 17 and the solution for the variational orbit through m 29. These powers have been chosen for another reason too. Double precision floating point arithmetic gives fifteen correct decimal places on our machine. Higher order terms in the series will not effect the fifteenth decimal digit of the answer when the series are evaluated at the given value for m. 2. The Intermediate Orbit

In Hill's lunar problem the Sun and. the Earth-Moon system rotate around eaizh other uniformly. The Sun is considered to be infinitely far away but the effect of its gravitation is still felt by the Earth-Moon system. If we introduce an earth centered coordinate system in which the positive x-axis always points in the direction of the 'fictitious' Sun then the equations for the position coordinates (x, y, z) of the Moon are x " - 2 n ' y ' - 3 n ' Z x = - t ~ r - 3x, y"

+ 2n'x'

Z tt "-}- H t 2 Z

=

-t~r-3y,

"-

~ [l,r- 3z.

In these equations r is the distance from the Earth to the Moon, n' is the mean motion of the Sun around the Earth-Moon system in an inertial coordinate system and/z is a constant which depends on the mass of the Earth-Moon system. Units are chosen in such a way that the gravitational constant is 1. Differentiation is with respect to time ' =

d/dt.

It is customary to introduce a new independent variable ~- by ~" = ( n -

n')(t-

to).

LITERAL SOLUTION FOR HILL'S LUNAR PROBLEM

281

n is the mean motion of the Moon around the Earth. The parameter to is unspecified and it can be used to determine different initial conditions on the same orbit. With m = n'/(n - n'), K = p / ( n - n ' ) 2 and . = d/dT the differential equations are 5~ - 2mj~ -

(2a)

3 m 2 x = - K x r - 3,

+ 2m2

(2b)

= -Kyr-3; + m 2 z = - K z r - 3.

.2

(2c)

We mention for later on that this system possesses the first integral j~2 _~. : 2

_~_ :~2 - - 3 m 2 x 2 _

2K/r - const.

For the following we can assume that K = 1 as it can be accomplished by a scaling of the form x ~ K1/ax, y ~ KI/ay, z -+ KX/3z. Hill starts the solution to the lunar problem by finding a 2rr-periodic orbit in the x - y plane which he calls the intermediate orbit. The solution ~ given as a power series in m with the coefficients only containing odd harmonic terms. To simplify the computations complex coordinates are introduced b~ u=

x + iy,

V-'-x--iy. ---- e i t .

With the operator d

d

D=~d~-

id-~.

the differential equations for the planar problem become DZu + 2 m D u

3

+Tm

D2v - 2 m D v + y3m

2'

(u + v) = U--1/2V--3/2~

2

(u + v) = u - 3 / 2 v -

1/2

(3)

.

It appears that Hill tried to avoid finding the expansion in m of the functions on the right hand side of the above equations and therefore used the so called homogeneous form of the differential equations which utilizes the existence of the first integral. The solution can then only be found up to a constant multiplicative factor. This factor is atso a power series in m and has to be found in a separate step. We propose to use (3) directly for finding u and its conjugate complex v. The expansions of the functions on the right hand side will be accomplished by an algorithm which is based on the following lemma. L E M M A 1: Given is r = 1 + ~

r,m". T h e series s = ~

s , m " which satisfies s = r '~

can be f o u n d recursively b y So = 1 a n d

s.=

(1 p=l

"'

n =

1,2,

....

282

Proof

D I E T E R S. S C H M I D T

The formula follows from the identity dr ds c~s dm = r drn"

We seek the solution to (3) in the form u =

Uo(r)

aktmk~2Z + 1

= k>_.O

0 ~ 2 Ill ~
The form of the solution will be justified by the fact that we can c o m p u t e the coefficients recursively with respect to k. The akz will be real rational numbers. F o r k = 0 the determining equation is aoo

=

--3 a oo

which has the solution aoo = 1. Assume now that the coefficients have been determined up to order k - 1. We will show how to find those of order k. Call ~:-lu = 1 + A and r = 1 + B. Then the coefficient of mkr 2l arising f r o m r can be written as dk~ - 1~.a kz _ ~ak,-z 3 where dkz is the corresponding coefficient from the series D = (1 + A)-1/2(1 + B) -3/2 -

1 + 89 + ~B.

Since this series starts with quadratic terms in m, dk~ will depend only on k n o w n aij with 0 ~< i < k. C o m p a r i n g coefficients of mk~ 2z+ 1 in (3) we obtain the determining relations ((2l + 1) 2 + 89 + ~3 a k . - z = Ok,, aak, + ( ( - - 2 l + 1)2 + 89 _ ~ = O k , - , for - k

(4) ~< 2l ~< k.

The right h a n d sides depend only on k n o w n quantities and are given by Ok, = dk, -- 2(2l + 1)ak_ ~., + 3 a k - 2., + 3 a k - 2, - , - ~ .

The restriction on I is due to the last term in Pk~ where we also used the convention that a~j = 0 for i < 0. F o r 1 = 0 (4) hag the solution ago = 89

and for l # 0 the solution is f

((21 ak l --

3

+

4/2(4/2 -

_,

1)

The m e t h o d which is described here for finding the coefficients akz follows closely the one described by Siegel and Moser (1970) in C h a p t e r 19. By majorizing the above series we can also present an easy convergence p r o o f similar to the one presented by Siegel. Unfortunately, the convergence radius which can be obtained by this method is very unsatisfactory as it is far below the one which can be obtained by more involved methods.

L I T E R A L S O L U T I O N FO R H I L L ' S L U N A R P R O B L E M

283

3. Hill's Differential Equation

The intermediate orbit of the last section is the reference orbit for the motion of the Moon in the theory of Hill and Brown. For m - 0 it is a circular Keplerian orbit. We will discuss the variation from the intermediate orbit which is due to the eccentricity and nonplanar motion of the Moon. In both cases a certain second order differential equation plays an important role. This equation is known as Hill's differential equation and has the form

(5)

ii + o(~)u = o.

0(r) is a Fourier cosine series containing only even harmonic terms. In solving Equation (5) Hill was led to the idea of an infinite determinant. When the computations are performed by machine another approach may be practical. We will use the Lie transformation for vector fields as it was developed by Henrard (1970) from the Lie transformation of Deprit (1969) for Hamiltonian functions. The method is summarized in the following lemma.

Given is the system of autonomous differential equations

L E M M A 2.

oo ~i i=0

i ' ~~ "

where e plays the role of a small parameter. Let the transJbrmation x = x(X, e) be the solution of the differential equation dx

de

= W(x, O

which satisfies the initial condition x ( X , O) = X. Assume that the generating Jhnction of this transformation has the form ,~=o i' W,+ ,(x). Then the transformed differential equation is

f(

"- ~

8i

,_-o i !

~'o(X)

and can be.found via a doubly indexed array of functions ~{ by i + k=O

The Lie derivative in the formula is given by

L ~ = e~ W ?x

~W~. ~x

In the lunar problem m plays the role of a small parameter and 0(r) in (5) has the

284

D I E T E R S. S C H M I D T

form 0(r) = ~ ' rn~0k(T) with 0o(~') = 1. Therefore we can treat (5) as a perturbation of a harmonic oscillator and solve it with the help of the Lie transformation in a way similar to the one described by Meyer and Schmidt (1978). Again it is convenient to introduce the complex variable ~: = e ~r and the operator D = - / ( d / d r ) . Hill's differential Equation (5) reads then D 2 u -- O(~)U = O,

(6)

where 0 has the special form Oo

0(~) = 1 + ~(~:)= 1 + ~

mk

k=l

~

Ak,~:2'

O~<21/l~
with Akl A k , - l" We transform (6) into a system and introduce complex coordinates by =

vl = u +

Du

v2 = u -

Ou

to obtain Dv,

= v, + l ( v ,

D v z = - v2 -

89

+ v2)O(~), + v2)O"(~),

(7)

DsC= ~. The third equation has been added to make the system autonomous and to stress the fact that ~ = e it. It should be noted that a real solution of (7) is given when Vl = ~2 and that for m - 0 the system represents three harmonic oscillators. In accordance with Lemma 2 set x = (vl, v2, ~:)r and find new coordinates X = (V~, I/2, ~:)r such that the transformed differential equations are as simple as possible. Because the third equation in (7) is already in normal form and because of the reality condition, the generating function of the transformation will be chosen in the form H' - -

W =~01

with

?Hi

y

~>_.o i!

Wi + I 9

At the different powers m i, i = 1,2, 3 , . . . , we will determine w~ such that the resulting differential equation is in normal form. All terms in the transformed equation can be eliminated except those which are in the kernel of Lwck ~ The first component of this kernel is determined by W-Vl

Ow Ow + v2 ~v~ ~v2

A typical term in w is cv~v[s I-

~

Ow -0. @

and it is in the kernel if and only if

c~+/3-7=0.

Since our system of differential equations is linear in vl and u2 only two combinations of exponents give terms which cannot be eliminated. The system (7) in normalized

LITERAL SOLUTION FOR HILL'S LUNAR PROBLEM

285

form is therefore

DVI = (1 + a)V~ + b~2Vz, DV2 = - b ~ - 2 V , - (1 + a)Z2, D ~ = ~,

(8)

where a and b are power series in m without constant terms. Another obvious transformation reduces this system to one with constant coefficients which is then used to find the solutions. When the term h = # a 2 - b 2 is different from zero, two linearly independent solutions of (8) are Vl = - b~:~ +x,

Vz = (a

1 +~

(9a)

Vl = - b ~ :~-x

Vz = (a + 1)~:-1-~

(9b)

-

~)~:-

and

Since ~: = e i~ the solutions are periodic if A is real. The Lie transformation also provides a m e t h o d of finomg the solution in the old vl, Vz variables. One transforms the functions vl, v2 or better yet at once the function u = 89 + Vz) by a m e t h o d which is similar to the one described in L e m m a 2. Instead of vector fields ~[ will represent functions ~o = 89 + v2), q~/o= 0, i = 1, 2, . . . are given and the transformed function ~ o~ (rni/i I ) ~ is found as in L e m m a 2 except that the Lie derivative is replaced by Lye

-

Ox W.

We obtain u = ~ mkck,~2'+1 V1 + ~ mkdk,~ 1t+1 V2, k,l

k l

where ck~ and dk~ are certain rational numbers. To obtain the solution to (5) we still have to replace V1 and V2 by (9a) or (9b) or a linear combination thereof. If we choose a linear combination which makes u real and which gives u = 0 for ~ = 0 we obtain the usual form of the solution co

u -- 2it=_~

Kt(~:zt+~_ ~:-2,-~).

(10)

Here we have set cr = h + 1 and K~ stands for a power series in m. In practice the normalization will be carried out only to a finite power in m say m N. In this case A can be found up to order N, but u only to order N - 1. The loss of one order is the fault of Equation (9) as the solutions Vl and Vz start with terms of order m. 4. Solution for the Inclination

As a direct application of the previous section we will find the first order a p p r o x i m a t i o n

286

DIETER S. SCHMIDT

for nonplanar orbits near the intermediate orbit. It requires that we find the solution to the variational equation of (2c). If ro = (UoVo)~/2 stands for the distance of the origin to the intermediate orbit, then the variational equation for the z-component is ;~'1 + ( m2 +

(11)

ro 3)zl = O.

This is Hill's differential equation with 0(r) = m 2 + r0-3 = 1 + ~ m k ~ l~
Mkzsezl

(12)

12/l~
where the given Mk~ satisfy Mk~ = M~,-z. Cowell (1896) used the method developed by Hill. He gave the solution to (11) as a power series with terms through m 6. The exponent o in (10) is denoted by g and Cowell computed terms up to m 8 for it. Nevertheless, the series expansions were not far enough to give the solution with the desired accuracy. Therefore CoweU and later on Brown used the numerical value (1) from the beginning to avoid the problem of the poor convergence of the series in m. With the help of a computer and with the method described in the previous section we can compute the series expansion as far as necessary. If we aim for the same accuracy as Brown (15 decimal places for g) we have to normalize (11) through order 16. The series expansion for g is given in Table I but the expansion of z l through m 15 is too long to be reproduced here. When the series are evaluated at the value (1) they agree with those of Brown's (1899). TABLE I Series e x p a n s i o n for g -- 1 + ~ gkm ~ k=l

k

gk

k

1

1.0

'9

2 3 4 5 6 7 8

0.75 -- 1.03125 --0.82031 0.02099 0.10445 0.58949 0.91020

25 60937 14973 61886 21676

5 96 94 52

10 11 12 13 14 15 16

gk

2.51936 7.36754 6.85832 --2.99983 -- 11.67490 --24.95418 -60.36517 - 109.96078

56212 66462 55862 75023 90071 85018 01179 50704

02 66 24 07

f

As the table shows the computation was performed with double precision floating point arithmetic. Although there would have been some advantages in doing the computations with rational arithmetic, we found that the integers which we had to store grew too large and exceeded the 15 decimal places, which is the maximum on our machine. The decision to stay with floating point arithmetic has another benefit. The method of the previous section works equally well when we start at Once with a special value

L I T E R A L S O L U T I O N FOR H I L L ' S L U N A R P R O B L E M

287

of m, that is we evaluate O(r) before we do the normalization. For reasons concerning programming it is beneficial to write (12) as O(r) = 1 + m ~ Ul~:2l,

where the N~ are now given numerically. The one m left in O(r) will play the role of the small parameter in the Lie transform. Thus again the solution is found as a series in m which is then evaluated at the given value of m. With this trick one computer program can either give the solution with m as a formal parameter or perform the computations with a special value of m. In the second case fewer terms in the series have to be computed for a desired accuracy, which is, of course, also the reason why Brown had chosen this method. 5. The Variational Equation w e return to the planar problem. The intermediate orbit in real form is x = Xo(

+ Vo('r))~

)

Y = Yo(') = 2-(Uo(,) - Vo(r)). Lt It satisfies (2a, b) which is an equation of the form "" - - 2m j? F~ y + 2m2 = Fy X

with F = r-1 + _}m2x2. The solution does not take into account the eccentricity e of the Moon's orbit. In a first approximation, one sets x - Xo .+ e x l , y = Yo + ey~ and determines xt and y~ from the variational equation ""

xl -

1 = Fx~xl + F~yYl,

Yl + 2ms

(13)

= F x , x l + Fy, y l .

These equations have a doubly periodic solution with periods 2rr and 2rr/e where c is unknown apriori. Before finding this solution it is common to determine e with the help of the differential equation for the normal displacement from the intermediate orbit (see Hill (1886); or any textbook i.e. Hagihara (1972)). I~et the tangent vector to the intermediate orbit be given by (Xo, Yo) -- V(cos 4', sin ~b). Then the tangential and normal displacement are related to the variation by a = xl cos 4~ + yl sin 4~, p = - x l sin4~ + Yl cos~. Due to the special nature of (13) the normal displacement satisfies HiU's differential equation +

o(.)p

= o

(14)

288

D I E T E R S. S C H M I D T

with

e(r) = i//V + 2(m + (b)2 + 2m z -- Fxx 89D (DZu~ D2 v~ Duo + Dv o ]

m

+ 89

+

Fyy

1 ~D2uo

4

D2uo

D2vo~

Duo

570)

Duo

DZvo~

2

+

2 m 2 mr-3"

If only isoenergetic variations are considered then the tangential displacement follows from d (-~)= dr

2(m + ~ ) ( - i f - ) .

(15)

The series in powers of m for c can 'be obtained from (14) but it converges very poorly. Hill (1894) found the series up to m 9 and through eleventh order in another parameter related to m. Hill's computations were checked and correctedby Bourne (1972). For practical purposes Hill and Brown used the numerical value (1) for rn to obtain more easily a numerical value for c which is correct to 15 decimal places. With it they returned to (13) to find the solution for x x and y~. T A B L E II Series expansion for c -- 1 +

~.~ c k m k k=l

k 1

2 3 4 5 6 7 8

9 10 11 12 13 14 15

ca 1.0 -O.75 -6.28125 - 1.84921 87500 00 E01 -5.45649 41406 25 E01 - 1.66666 30045 57 E02 - 5 . 6 3 7 8 1 80101 18 E02 - 2.13428 24002 02 E03 - 8 . 7 9 4 0 1 47141 78 E03 -3.79648 58198 74 E04 - 1.67291 46049 10 E05 - 7 . 4 4 5 7 9 18137 01 E05 - 3 . 3 4 0 0 3 70913 78 E06 - 1.51065 90916 51 E07 - 6 . 8 9 1 6 2 86361 17 E07

k

ck

16 17

--3.17023 43642 E08 - 1.46936 73423 E09 - 6 . 8 5 5 2 4 16886 E09 - 3 . 2 1 6 4 9 43122 E l 0 - 1.51668 00033 E l l - 7 . 1 8 3 0 4 76725 E l l - 3 . 4 1 5 3 1 87281 E l 2 - 1 . 6 2 9 6 7 37673 E l 3 - 7 . 8 0 1 5 3 17158 E l 3 - 3.74583 61071 E l 4 - 1 . 8 0 3 4 3 72535 E l 5 - 8 . 7 0 4 4 6 81480 E l 5 - 4 . 2 1 1 0 3 11744 E l 6 - 2 . 0 4 1 5 7 34453 E l 7 - 9 . 9 1 7 5 6 99800 E l 7

18

19 20 21 22 23 24 25 26 27 28 29 30

p w

If we want the same accuracy from a literal expansion for c we have to compute terms through m 3~ The terms of the series are given in Table II. For the computation of 0(~-) we used our solution Uo(T) for the intermediate orbit which contains terms through order 18 in m. Although 0(~-) is expanded through order 30 and the Lie transformation is carried to the same order, the use of a truncated series Uo(~-) does not impair the desired accuracy. As mentioned in Section 3 our procedure provides without much additional effort the series expansion for p. We use it to determine or/V from (15). This requires only a

LITERAL SOLUTION FOR HILL'S LUNAR PROBLEM

289

TABLE III Corrections to Brown's solution of the variational equations. Brown (1899, page 94).

Corrections to Brown's value

New value

f

e_ 1 -0.14889 75285 5 e_2 -0.00005 20895 6

+ 0.00000 00011 - 0.00000 00032

straightforward integration after the known together.lThe only series which has not been which has to be found from V - ~ = ( - D u o To obtain x l and Yx only one additional one notices that Ul

= xl

series on the right of (15) are multiplied used explicitly earlier is the one for V-1 D v o ) - 1 / 2 and the help of Lemma 1. multiplication of series is required after

/ + iyx = ((7 + ip)e i~ = ( i \

(7

V

Again the series are too long to be reproduced here since at order m k there are 2k + 2 terms. But we did evaluate the series at m and compared it with Brown's result. He gives the solution for u l in the form ul~ -1 = ae ~ ( e , ~ :2'+c + e ~ 2 ' - c ) , i

where the parameters a and e are chosen in such a way that to - to = 1. The coefficients are listed with 10 decimal places. In contrast to the solution for the inclination we found here some differences between ours and Brown's solution. The two coefficients with the largest discrepancies are listed in Table III. We double checked our result by redoing the computations with the special value of m by the method which was described at the end of Section 4. References Bourne, S. R.: 1972 Celest. Mec. 6, 167-186. Brown, E. W. : 1899, Mem. Roy. Astron. Soc. 53, 39-116, 163-202. Brown, E. W. : 1904, Mem. Roy. Astron. Soc. 54, 1-64. Brown, E. W.: 1905, Mem. Roy. Astron. Soc. 57, 1-145. Brown, E. W.: 1908, Mem. Roy. Astron. Soc. 59, 1-103. Cowell, P. H.: 1896, Am. J. Math. 18, 99-127. Deprit, A. : 1969, Celest. Mech. 1, 12-30. Hagihara, Y. : 1972, Celestial Mechanics, Vol. 2, Part 2, MIT Press. Henrard, J.: 1970, Celest. Mech. 3, 107-120. Hill, G. W.: 1878, Am. J. Math. 1, 5-26, 129-147, 245-260. Hill, G. W. : 1886, Acta. Math. 8, 1-26. Hill, G. W. : 1894, Ann. Math. 9, 31-41. Knuth, D. E.: 1968, The Art of Computer Programming, Vol. 3, Addison-Wesley. Meyer, K. R. and Schmidt, D. S.: 1977, Funkcialaj Ekvacioj 20, 171-192. Rom, A.: 1969, Celest. Mech. 1, 301-319. Rom, A.: 1971, Celest. Mech. 3, 331-345. Siegel, C. L. and Moser, J. K." 1970, Lectures on Celestial Mechanics, Springer.