The Gravitational Ellipse

JOURNAL OF MATHEMATICAL PHYSICS 50, 022901 共2009兲

The gravitational ellipse Maurizio M. D’Eliseoa兲 Osservatorio S. Elmo, Via A. Caccavello 22, 80129 Napoli, Italy 共Received 27 July 2008; accepted 12 January 2009; published online 19 February 2009兲

The elliptical orbit of the classical gravitational two-body problem can be determined by studying the free oscillations about a circular motion or the small motions around a fixed point in a rotating reference frame. In this last schematization we approximate the differential equation of motion by a succession of simple equations we solve iteratively, obtaining a piecemeal determination of the position vector r formally expressed in terms of Laurent polynomials, from which we quickly deduce the explicit time-dependent expressions in the form of complex trigonometric polynomials. This approach can also be used in the presence of perturbing forces and, by way of illustration, we study the effects of a small linear repulsive force on the elliptical orbit. © 2009 American Institute of Physics. 关DOI: 10.1063/1.3078419兴

I. INTRODUCTION

In the classical gravitational two-body problem we can easily derive from Newton’s equation of motion, by means of the area and Laplace integrals, the elliptical shape of the orbit in function of the polar angle ␪.1,2 In the practical applications we look at the position of the orbiting body in function of time, but the time dependence of this variable cannot be expressed by closed-form functions. One way out is to obtain, starting from the area integral, by means of repeated approximate integrations, the function ␪共t兲 in the form of a cosine series whose coefficients are, in turn, power series in the eccentricity. Then, following the route t → ␪共t兲 → r共␪兲, we can at last find the position vector r for any time to arbitrary accuracy. An alternate approach, which we will follow in this paper, is to obtain the series development of r共t兲 by a direct integration of the equation of motion, opportunely modified with a transformation of the force term. Central to our analysis is the framework of complex variables because of the following. • From the form of the sought solution, we are led to consider a uniformly rotating reference frame in the plane of motion, and we can formally represent the gravitational force as a sum of homogeneous polynomials in two complex conjugate variables. • We can substitute the two scalar equations of this dynamical system with two degrees of freedom, by a single complex equation, containing on the left-hand side a linear differential operator with constant coefficients, and on the right-hand side nonlinear source terms having small parameters, namely the powers of the eccentricity, as factors. Such a form allows the application of the method of iteration, and the path r共t兲 is then determined solving a succession of elementary linear equations by a well defined procedure. • Knowing r共t兲, expressed in complex form, to as high a degree of precision as we wish, we can easily deduce, by routine methods of complex algebra and calculus, the related time developments of distances, angles, and of all other orbital functions that are of interest in astrodynamics. • We may employ a simple notation for the complex exponential, which will ease all computations, either done by hand or by computer.

a兲

Electronic mail: [email protected].

0022-2488/2009/50共2兲/022901/10/$25.00

50, 022901-1

© 2009 American Institute of Physics

Downloaded 26 Feb 2009 to 192.135.20.186. Redistribution subject to AIP license or copyright; see http://jmp.aip.org/jmp/copyright.jsp

022901-2

J. Math. Phys. 50, 022901 共2009兲

Maurizio M. D’Eliseo

y y’

B

x’

as

r a nt

A

x

FIG. 1. Geometry of motion.

II. THE EQUATION OF MOTION

Given the two point masses 共A , mA兲 and 共B , mB兲 interacting gravitationally, the equation of motion of B in an A-centered reference system in the orbital plane is1,2 r¨ = − ␮r−3r,

共1兲

where r = x + iy = rei␪ is the complex coordinate of B, r2 = 兩r兩2 = rrⴱ, ␮ = G共mA + mB兲, G is the gravitational constant. This equation admits the circular solution r = aeint of angular frequency n = 冑␮a−3 and satisfying the initial condition r共0兲 = a. In order to find the general periodical solution of 共1兲 after a transformation of the awkward r−3 term, we first put the equation into a more convenient form by replacing the independent variable time t with the scaled imaginary time ␶ ⬅ int. Accordingly, denoting by a prime ⬘ the differentiation with respect to ␶, Eq. 共1兲 turns into r⬙ = a3r−3r.

共2兲 ␶

It is helpful to write the circular solution in the form r = a␭, where ␭ ⬅ e = e = cos nt + i sin nt. For any integer k the conjugate of ␭k satisfies the identity 共␭k兲ⴱ ⬅ ␭−k, while the derivative of order j amounts to a multiplication by k j, so we can regard the differentiation sign ⬘ acting on integer powers of ␭ as the logarithmic derivative d / d ln ␭ or the Euler operator D ⬅ ␭ · d / d␭. This notation hides the presence, in the formulas, of both the time and the imaginary unit, and we can work with simple algebraic objects, the Laurent polynomials in ␭ on the unit circle 兩␭兩 = 1, which exploit the identity between conjugate and reciprocal of the complex exponential. After computing the polynomials that are solutions of our problem we will replace ␭ by eint. We study now the near-circular solutions of 共2兲 by assuming r = a共1 + s兲␭, where s = s共␭兲 is a complex function3 whose range is the open unit disk d共0 ; 1兲. We take s complex in order to allow for motions along arbitrary directions in the plane. Note that being s the composite function s共␭共␶共t兲兲兲, since ␭ = 1 when t = 0, the initial value of s will be s共1兲. Insertion of the assumed expression of r into 共2兲 and cancellation of the common factor a␭ yield D2s + 2Ds + s + 1 = 共1 + s兲−1/2共1 + sⴱ兲−3/2 .

int

共3兲

This equation describes the small motions of B around the point lying on the rotating radius vector of the circular orbit scaled at unit distance to the origin. Besides, in 共3兲, from left to right there are terms corresponding to the forces acting: inertial, Coriolis, centrifugal, and gravitational. All this highlights that, from a physical point of view, the equation is related to a frame x⬘y ⬘ uniformly rotating with angular velocity n 共see Fig. 1兲. In this reference system the circular solution is s = 0. It is interesting to notice that the rotating system naturally stems from our treatment of the near-circular motion in the framework of complex variables. Knowing s, the resulting motion in the nonrotating inertial frame is found by computing r = a共1 + s兲␭, from which we can immediately write out the position vector r共t兲 = aeint + as共t兲eint.

Downloaded 26 Feb 2009 to 192.135.20.186. Redistribution subject to AIP license or copyright; see http://jmp.aip.org/jmp/copyright.jsp

022901-3

J. Math. Phys. 50, 022901 共2009兲

The gravitational ellipse

We get rid of the irrational factors in the right-hand side of 共3兲 by developing them, in turn, in formal power series in s and sⴱ, performing their Cauchy product, and writing on the left-hand side the two resulting linear terms,

共4兲 It is desirable to have an equation for s only because, knowing s, we also know its conjugate sⴱ. This can be done by the formal device of introducing the star ⴱ conjugation operator, with the property that ⴱX = Xⴱ, where X is any complex expression. Then, by defining the linear differential operator Lˆ ⬅ D2 + 2D + 23 + 23 ⴱ ,

共5兲

Lˆs = S.

共6兲

we cast 共4兲 in the form

ⴱ

The nonlinear source function S = S共s , s 兲, as shown in 共4兲, can be regarded as an infinite sum of homogeneous polynomials in s , sⴱ of increasing degree, which in practice we truncate somewhere, so that 共6兲 will have the form K

K

␯

Lˆs = SK ⬅ 兺 S共␯兲 ⬅ 兺 兺 ␴␯−j,js␯−jsⴱj , ␯=2

␯=2 j=0

共7兲

where the coefficients ␴␯−j,j are rational numbers that can be expressed by means of double and/or simple factorials,

␴␯−j,j = 共− 1兲␯

关2共␯ − j兲 − 1兴 ! ! 共2j + 1兲 ! ! 关2共␯ − j兲兴 ! 共2j + 1兲! = 共− 1兲␯ . 2␯共␯ − j兲 ! j! 关2␯共␯ − j兲 ! j!兴2

共8兲

We recall that the double factorial j ! ! is recursively defined: j ! ! = 1, if j = −1 , 0 , 1; j ! ! = j共j − 2兲 ! !, if j ⱖ 2, while 共2j − 1兲 ! ! = 共2j兲 ! / 2 j j ! , 共2j + 1兲 ! ! = 共2j + 1兲 ! / 2 j j!. In 共7兲, the homogeneous polynomial S共␯兲 represents the block gathering the ␯ + 1 terms of the source of the same degree ␯, while the upper index K corresponds to the overall approximation degree, which we fix in advance. III. THE MAKING OF THE ELLIPSE

The integration of Eq. 共7兲 can be effected by means of a straightforward iterative scheme, which can be described as consisting in the construction and in the solution of a sequence of equations ever more approximating the equation of motion 共6兲. If S j stands for S共s j , sⴱj 兲, we put Lˆs1 = 0, j , Lˆs j = S j−1

j = 2, . . . ,K.

共9兲 共10兲

At each step we build a particular right-hand source, and because of the nature of the former solutions, the corresponding equation becomes a linear differential equation whose solution is the sum of the general integral s1 of the homogeneous equation 共9兲 and of the particular integrals due to the presence of the source terms. In practice, by the linearity of Eq. 共10兲, we can decouple them, and to find s j we will identify every time in both sides only the terms of the same order 共in the eccentricity, as we will shortly see兲. This can be easily done for high j’s by computer taking the difference of the power series

Downloaded 26 Feb 2009 to 192.135.20.186. Redistribution subject to AIP license or copyright; see http://jmp.aip.org/jmp/copyright.jsp

022901-4

J. Math. Phys. 50, 022901 共2009兲

Maurizio M. D’Eliseo

developments in the eccentricity of the source function to the orders j and j − 1, respectively. The rationale is that the solutions concerning the terms of order less than j have supposedly already been found, while those related to the terms of order higher than j would be incomplete because of the intrinsic nature of the approximation process and the formal structure of the source terms. Therefore it suffices to look for the particular integral of the equation relative to the unknown function s共j兲 ⬅ s j − s j−1, of order j in the eccentricity, and to substitute 共10兲 with the reduced scheme, j 兴, Lˆs共j兲 = 关S j−1

j = 2, . . . ,K,

共11兲

where we agree that the square brackets enclose all and only the terms of order j of the forcing function. Let us start from 共9兲 and, in full generality, assume as solution the binomial s1 = a␭k + b␭−k, where the positive integer k and the complex constants a , b have to be determined. After developing the expression Lˆ共a␭k + b␭−k兲, equating separately to zero the resulting coefficients of ␭k , ␭−k, and taking the conjugate of the second equation, we shall obtain a homogeneous linear system in the unknowns a , bⴱ, which may be written as the matrix equation

冉

k2 + 2k + 3/2 3/2

3/2 2

k − 2k + 3/2

冊冉 冊 冉冊 ·

0 a , ⴱ = 0 b

共12兲

which becomes singular, admitting so nontrivial infinitely many solutions, provided k = 1, since its determinant is k4 − k2,4 and accordingly fixing to n the frequency of the homogeneous solution. It follows the relation 3a = −bⴱ from which, if we let a = 共1 / 2兲eⴱ, we obtain b = −共3 / 2兲e. In hindsight, we select e = eei␻, the eccentricity vector, where 0 ⬍ e ⬍ 1 is the 共scalar兲 eccentricity and ␻ is the argument of pericenter.1,2 Then the solution of 共9兲 is s1 = 21 eⴱ␭ − 23 e␭−1 ,

共13兲

r1 = a共1 + s1兲␭ = − 23 ae + a␭ + 21 aeⴱ␭2 = − 23 ae + aeint + 21 aeⴱe2 int .

共14兲

while in the inertial frame

We may verify at once that 共14兲 agrees with the expression, in function of time, of the elliptical orbit to order e obtained with the usual methods in complex form.1,2 Starting from r = aei␪共1 − e2兲 / 关1 + e cos共␪ − ␻兲兴 and from the area integral ᐉ = 冑␮a共1 − e2兲 = na2冑共1 − e2兲, we have, to order e, r = rei␪ ⬇ a关1 − Re共eⴱei␪兲兴ei␪ = − 21 ae + aei␪ − 21 aeⴱe2i␪ ,

共15兲

ᐉ na2 ␪˙ = 2 ⬇ 2 ⬇ n + 2en cos共␪ − ␻兲, r r

共16兲

␪ ⬇ nt + 2e sin共nt − ␻兲 = nt − i共eⴱeint − ee−int兲.

共17兲

From 共17兲 it is easily seen that for a first-order accuracy we can assume ei␪ ⬇ − e + eint + eⴱe2 int,

e2i␪ ⬇ e2 int ,

共18兲

and, after substituting the exponentials of 共18兲 in 共15兲, the expression 共14兲 of r1 is retrieved. Before go on, we illustrate the key points of the procedure to be followed. • To simplify the computations, in s1 we will put e , eⴱ = e 共␻ = 0兲, so that, when t = 0, we shall have s1共1兲 = −e, and B will be at pericenter on the real axis, since r1共0兲 = a共1 − e兲. Analogously, when t = ␲ / n, ␭ = −1, s1共−1兲 = e, r1共␲ / n兲 = a共1 + e兲, and B will be at apocenter. This choice of the initial conditions for s1 shall reverberate on the successive solutions s共j兲 , j ⬎ 1,

Downloaded 26 Feb 2009 to 192.135.20.186. Redistribution subject to AIP license or copyright; see http://jmp.aip.org/jmp/copyright.jsp

022901-5

J. Math. Phys. 50, 022901 共2009兲

The gravitational ellipse

and the consequent conditions s共j兲共⫾1兲 = 0 will be automatically satisfied for all even j ⱖ 2, while, in order to assure the validity of the same conditions for the odd-j solutions, it will be required we fix, in them, to zero the value of the coefficients of ␭−1. We will restore the full form of the vectors e , eⴱ at the end since, given the initial couplings e␭−1 ↔ e␭−1 , eⴱ␭ ↔ e␭, in the results we shall find terms of the form e␮␭⫾␯ 共␮ , ␯ having the same parity and ␮ ⱖ ␯兲, which we will translate back as e ␮␭ −␯ → e 共␮−␯兲e ␯␭ −␯,

e ␮␭ ␯ → e 共␮−␯兲e ⴱ␯␭ ␯ , i␯␻

共19兲

−i␯␻

and so we shall recover the phase factors e , e , characterizing a generic pericenter position, whence an arbitrary orientation of the orbit in the plane. This operation should be effected on sK before returning to the inertial frame. • Doing the iterations outlined in 共11兲, we shall encounter equations of the form Lˆs共j兲 = 兺 共pi␭i + qi␭−i兲,

i = j, j − 2, . . . ,

i ⱖ 0,

共20兲

i

where 2 ⱕ j ⱕ K and i runs through non-negative integers with the same parity of j. The coefficients pi , qi are rational monomials in e j, and it happens that always p1 is a multiple of q1, and precisely of the form p1 = 3q1. It is worth noting that in 共20兲 the forcing function for ␭ = 1 verifies the identity 兺i共pi + qi兲 = 共j + 1兲e j. We have here a manifestation of the inversesquare force at work: when ␭ = 1, at pericenter, it follows that a2r−2 = 共1 − e兲−2 = 兺 j共j + 1兲e j. An analogous relation with alternating signs is found for ␭ = −1 at apocenter. This enables one to do intermediate checks during the computations. To solve the Eq. 共20兲, we will apply the method of undetermined coefficients,4 by assuming as solution the sum of a function of the same type of the known source and of all its independent derivatives, so that it suffices to deal with the following three types of equations: Lˆa0 = p0 ,

共21兲

Lˆ共a1␭ + b1␭−1兲 = p1␭ + q1␭−1, Lˆ共ak␭k + bk␭−k兲 = pk␭k + qk␭−k,

p1 = 3q1 ,

共22兲

k ⬎ 1,

共23兲

where the constants a’s and b’s have to be found. From 共21兲 we immediately find 3a0 = p0 → a0 = 31 p0 .

共24兲

In order to satisfy 共22兲, because of the relation existing between the coefficients of the source function, either a1 or b1 may be chosen arbitrarily, but the two constants must satisfy the relation 3a1 + b1 = 2q1. We try b1 = 0, as to have Lˆ共a1␭兲 = 3q1␭ + q1␭−1 ⇒ 3a1 = 2q1 ⇒ a1 = 32 q1 .

共25兲

It turns out that this choice ensures the fulfillment of the initial conditions s j共⫾1兲 = 0 for all odd j ⱖ 3. Last, we consider 共23兲. From 共12兲, with a = ak , bⴱ = bk, it follows that the real coefficients ak , bk , pk , qk satisfy to the matrix equation

冉

k2 + 2k + 3/2 3/2

3/2 2

k − 2k + 3/2

冊冉 冊 冉 冊 ·

ak pk = , bk qk

which we can write in abridged form

Downloaded 26 Feb 2009 to 192.135.20.186. Redistribution subject to AIP license or copyright; see http://jmp.aip.org/jmp/copyright.jsp

022901-6

J. Math. Phys. 50, 022901 共2009兲

Maurizio M. D’Eliseo

Lk · 共ak,bk兲T = 共pk,qk兲T ,

共26兲

where the transposition symbol T transforms a row vector into a column vector. Since det共Lk兲 ⫽ 0, the solution is found multiplying on the left both sides of 共26兲 by the inverse matrix L−1 k , so we have T 共ak,bk兲T = L−1 k 共pk,qk兲 ,

共27兲

and it follows that the coefficients ak , bk are given by

ak =

2共k2 − 2k兲pk + 3共pk − qk兲 , 2共k4 − k2兲

bk =

2共k2 + 2k兲qk − 3共pk − qk兲 . 2共k4 − k2兲

共28兲

• At the end of the computations, we shall obtain sK, the sum of the partial solutions of increasing orders, expressed as a Laurent polynomial, which will satisfy the pericenter/ apocenter conditions, in the form K

sK = s1 + 兺 s共j兲 = j=2

+K

兺

f i␭ i,

sK共⫾1兲 = ⫿ e,

共29兲

i=−K

where the coefficients f i are, in turn, rational polynomials in e, either odd or even according to the parity of i. In conclusion, the solution sK of 共7兲 is of order eK, has in the real time the overall period T = 2␲ / n, and in the inertial frame the vector rK = a共1 + sK兲␭, likewise algebraically expressed as a Laurent polynomial in ␭, will have the same period. Working through the scheme 共11兲 up to order e3 we successively get the equations Lˆs共2兲 = − 23 e2 +

Lˆs共3兲 = −

27 3 16 e ␭

−

15 2 2 4e ␭

9 3 −1 16 e ␭

+

+ 43 e2␭−2 ,

89 3 3 16 e ␭

+

11 3 −3 16 e ␭ .

Following the rules just described, we find s3 = s1 + s共2兲 + s共3兲 = 21 e␭ − 23 e␭−1 − 21 e2 + 83 e2␭2 + 81 e2␭−2 − 83 e3␭ + 31 e3␭3 +

1 3 −3 24 e ␭ .

共30兲

As last step, after applying the translation keys 共19兲, we compute r3 = a共1 + s3兲␭, and, turning back to the exponential notation, we get the complex trigonometric polynomial

冉

冊

冉

冊

1 1 3 1 3 1 3 1 r3 = − e + 1 − e2 eint + e2e−int + eⴱ − e2eⴱ e2 int + e3e−2 int + eⴱ2e3 int + eⴱ3e4 int . 2 2 8 2 8 24 8 3 a 共31兲 Whatever the value of K, we can verify that the orbit described by rK is the approximate expression of an ellipse. Let us recall the ellipse’s standard definition: it is a locus of points in a plane such that the sum of the distances to two fixed points 共the foci兲 is a constant 共given by twice the semimajor axis a兲. In our case, the first focus is located in the origin, where body A rests, while the second and empty focus resides in the point −2ae, and we can verify that the sum of the two focal distances, to order eK, is 兩rK兩 + 兩rK + 2ae兩 = 2a, as it should be for a true ellipse. Also, if ␥K is the closed path described by rK in one period, complex analysis5 tells us that the enclosed area is AK = 共1 / 2i兲养␥KrKⴱ drK, and we actually find

Downloaded 26 Feb 2009 to 192.135.20.186. Redistribution subject to AIP license or copyright; see http://jmp.aip.org/jmp/copyright.jsp

022901-7

J. Math. Phys. 50, 022901 共2009兲

The gravitational ellipse

AK =

1 2i

冕

T

0

冉

冊

1 1 rKⴱ 共t兲r˙ K共t兲dt = ␲a2 1 − e2 − e4 − ¯ , 2 8

共32兲

which agrees with the ellipse area ␲ab ⬅ ␲a2冑1 − e2 to the order eK. We have thus compelling reasons to guess and state that the curve described by rK for K → ⬁ is an ellipse, even if we had not known in advance by other means the true nature of the orbit. The expression of rK共t兲 is a parametric representation of the ellipse traveled by respecting the area’s law. From it we can obtain, by means of the usual operations on complex expressions and some power series developments, all information, to order eK, about the orbiting body, as distances, angles, velocities, expressed in any functional form, in function of time. We get so all the expansions of the elliptical motion of celestial mechanics.6 As an example, we deduce the expressions of the rectangular and polar coordinates. At once we find x = Re共r兲 = Re关a共1 + s兲␭兴 = a cos nt + a Re关s共t兲eint兴,

共33兲

y = Im共r兲 = Im关a共1 + s兲␭兴 = a sin nt + a Im关s共t兲eint兴,

共34兲

besides, since r = rei␪ = a关1 + s共t兲兴eint , r2 = a2共1 + s + sⴱ + ssⴱ兲, from the series developments in s , sⴱ we get r ⬇ a + a Re共s − 41 s2 + 41 ssⴱ + 81 s3 − 81 s2sⴱ − ¯兲 ,

共35兲

␪ = − i ln共ei␪兲 = − i ln共r/r兲 = − i ln关eint共1 + s兲/兩1 + s兩兴 ⬇ nt + Im共s − 21 s2 + 31 s3 − ¯兲 .

共36兲

IV. A LINEAR REPULSIVE FORCE

If the conditions of the problem are slightly modified by the presence of some perturbation term, the method under consideration can be applied without any alteration. Let a planet P⬘ of mass m⬘ be moving along the outer circular orbit r⬘ = a⬘ei␪⬘ ⬅ a⬘ein⬘t, whose angular frequency n⬘ is incommensurable with that n of planet B, whose orbit lies in the same plane. We will consider the direct force F of P⬘ acting on a generic fixed point r = rei␪ of the orbit of the planet B. We have also an indirect force term due to the influence of P⬘ on B through the perturbation of the central body A. This force is ␮⬘r⬘ / r⬘3, but since its average value reduces, by the formulas that shortly shall follow, to the computation of the integral 兰20␲ei␪⬘d␪⬘ = 0, it does not concern us here. Let now

␣ = r/a⬘ Ⰶ 1, so that d␾ = d␪⬘. Setting ␮⬘ ⬅ Gm⬘ ,

⌬ = 1 + ␣2 − 2␣ cos ␾,

␾ = ␪⬘ − ␪ ,

⑀ ⬅ ␮⬘ / a⬘3, we have

F = ␮⬘

共a⬘ei␾ei␪ − r兲 r⬘ − r = ⑀ . ⌬3/2 兩r⬘ − r兩3

共37兲

The average value of F over the period T⬘ of P⬘ can be expressed in the form 具F典 =

1 T⬘

冕

T⬘

0

F dt =

1 2A⬘

冕

A⬘

F dA⬘ ,

共38兲

0

where we have used the area integral and the circumstance that the area’s constant is twice the areal velocity. In this expression we have the quotient of two areas: the area dA⬘ / 2 = a⬘2d␪⬘ / 2 of the elementary sector with the vertex in the central body and the total area A⬘ = ␲a⬘2 of the orbit of P⬘. If we set

Downloaded 26 Feb 2009 to 192.135.20.186. Redistribution subject to AIP license or copyright; see http://jmp.aip.org/jmp/copyright.jsp

022901-8

J. Math. Phys. 50, 022901 共2009兲

Maurizio M. D’Eliseo

y y’ B as

r

r’

a

nt

A

P’

x’

n’t

x

FIG. 2. External planet perturbation.

d␮⬘ ⬅

␮⬘dA⬘ ␮⬘a⬘2d␪⬘ ␮⬘d␪⬘ = = , 2 ␲ a ⬘2 2␲ 2A⬘

共39兲

we may write 共38兲 in the form 具F典 =

冖

r⬘ − r d␮⬘ , 兩r⬘ − r兩3

共40兲

where the loop integral is taken along the orbit of P⬘ in the direction of motion. We suppose that the two orbits do not intersect each other. This result demonstrates that 具F典 is equivalent to the force exerted, on the point r, by a thin homogeneous material wire ring of total mass m⬘, which can be considered as an ideal materialization of the circular orbit of P⬘. The averaging process is meaningful since, by the supposed incommensurability of the two mean motions n , n⬘, at any time the two positions r, r⬘ are stochastically independent. If the motion of P⬘ is elliptical and noncoplanar to that of B, then r , r⬘ are ordinary vectors, and we shall have, in formula 共39兲, dA⬘ = r⬘2d␪ / 2 , A⬘ = ␲a⬘b⬘ = ␲a⬘2冑共1 − e⬘2兲. The distribution of the mass on the now elliptical ring is no more uniform: the mass of any element of the ring is proportional to the time this element is traversed by the perturbing planet. This theorema hoc elegans was the starting point of an important work by Gauss on the secular planetary perturbations.7 From 共38兲 we now compute the average force acting on the point r, and to do this it must be considered that 1 / T⬘ = n⬘ / 2␲, while, from the area integral for the circular orbit of P⬘ given by a⬘2␪˙ ⬘ = n⬘a⬘2, it follows that dt = d␪⬘ / n⬘ = d␾ / n⬘. Then 具F典 =

1 2␲

冕

2␲

Fd␾ =

0

⑀ 2␲

冕冉 2␲

0

冊

r i␾ d␾ ⑀a⬘ r a⬘e − r 3/2 ⬇ r ⌬ 2␲ r

冕

2␲

0

共ei␾ − ␣兲共1 + 3␣ cos ␾ + ¯兲d␾

1 = ⑀r + O共␣2兲. 2

共41兲

Thus, to order ␣, an external planet moving on a circular orbit 共as well as a circular homogeneous material ring兲 exerts, approximately, a linear repulsive force, directed away from the center, on a particle located somewhere near the center 共see Fig. 2兲 and traveling along a coplanar orbit. To insert this force in our equation of motion 共6兲, we substitute, in 共41兲, r by a共1 + s兲␭, divide by −n2a␭, and we get, in the rotating frame, Lˆs = G共s,sⴱ兲 ⬅ − 21 w2 − 21 w2s + S共s,sⴱ兲,

共42兲

where w ⬅ 冑⑀ / n Ⰶ 1. We could integrate this equation iteratively, starting from the circular orbit 共s0 = 0兲, but the first and only equation we need to consider here is

Downloaded 26 Feb 2009 to 192.135.20.186. Redistribution subject to AIP license or copyright; see http://jmp.aip.org/jmp/copyright.jsp

022901-9

J. Math. Phys. 50, 022901 共2009兲

The gravitational ellipse

Lˆs1 = − 21 w2 → s1 = − 61 w2

共43兲

because it is evident that further iterations will produce a succession of constant terms of increasing powers of w, leading to negligible additional effects. Equation 共43兲 describes a slight shrinking of the unperturbed circular orbit of B, as we have r = a共1 − w2 / 6兲␭. We study now the small motions ␦s about this modified circular orbit, and to do this we consider the variation equation deduced from 共42兲,

␦共Lˆs − G兲 = Lˆ␦s − 共⳵sG兲s1␦s − 共⳵sⴱG兲s1␦sⴱ = 0,

共44兲

where we set, to the lowest order in s,

冉 冊 ⳵G ⳵s

1 3 3 3 ⬇ − w2 − s1 − sⴱ1 − ¯ = − w2 , 2 4 4 4 s1

冉 冊 ⳵G ⳵ sⴱ

3 15 3 ⬇ − s1 − sⴱ1 + ¯ = − w2 . 4 4 4 s1

If we put Lˆ⬘ ⬅ D2 + 2D + W + W ⴱ ,

W ⬅ 23 + 43 w2 ,

共45兲

Equation 共44兲 may be written in the form Lˆ⬘␦s = 0.

共46兲

We observe that, in the limit w → 0, the operator Lˆ⬘ goes to Lˆ, and Eq. 共46兲 matches the homogeneous Eq. 共9兲 of the unperturbed motion. Knowing that the solutions of well behaved differential equations depend continuously on the coefficients, we assume the following form of the solution: 1 ␬ ␦s = heⴱ␭c + he␭−c , 2

2

共47兲

where the real constants h, c ⬇ 1, and ␬ ⬇ −3 for w → 0. The constant h represents a slight correction to the value of the scalar e, while c implies a uniform motion of the perihelion, since we may write e␭−c = e␭1−c␭−1 = eei关␻+共1−c兲nt兴␭−1 ⬅ e+␭−1 ,

共48兲

and the perihelion advances or recedes according as c ⭴ 1. Plugging 共47兲 into 共46兲, we obtain the nonlinear system of equations in c , ␬, c2 + 2c + 共1 + ␬兲W = 0,

共49兲

␬共c2 − 2c + W兲 + W = 0,

共50兲

which can be managed by putting, in 共50兲, c = 1, solving for ␬

␬=

W ⬇ − 3 + 3w2 , 1−W

共51兲

and inserting this value of ␬ in 共49兲 so as to obtain in the second approximation,

Downloaded 26 Feb 2009 to 192.135.20.186. Redistribution subject to AIP license or copyright; see http://jmp.aip.org/jmp/copyright.jsp

022901-10

J. Math. Phys. 50, 022901 共2009兲

Maurizio M. D’Eliseo

3 3 ⑀ c = − 1 + 冑1 − 共1 + ␬兲W ⬇ 1 − w2 = 1 − . 4 4 n2

共52兲

The determination of the value of the last constant h depends from the meaning we assign to the eccentricity e in the perturbed motion. We leave aside the discussion on this matter and simply point out that, by the supposed smallness of ⑀, we can safely assume h = 1. Then the solution we have found is, in the inertial frame, r1 = 共1 + s1 + ␦s兲␭ = 共1 − 61 w2兲␭ + 21 aeⴱ␭1+c − 23 ae␭1−c = − 23 ae+ + a共1 − 61 w2兲eint + 21 ae+ⴱe2 int , 共53兲 to be compared with Eq. 共14兲 of the unperturbed motion. Here the perihelion moves at the constant rate 共1 − c兲n while, considering the relation n = 冑␮a−3 which is valid for the unperturbed circular orbit, one can say that, in the actual orbit, for a given mean motion n the mean distance is smaller than it would be without disturbance. Let us do now a quick application to the planet Mercury perturbed by Saturn, neglecting the eccentricity of the latter and the mutual inclination of the orbital planes. We consider first, using 共48兲, the variation in the argument of Mercury’s perihelion after the completion of a single orbit, ⌬␻ = 共1 − c兲n⌬t, where ⌬t = T = 2␲ / n and c has the form 共52兲 with

冉

3 ⬘ amer 1 0.387u.a. ⑀ Gmsat = 3 = 2 3497.9 9.582u.a. n Gmsun asat ⬘

冊

3

= 1.883 ⫻ 10−8 .

共54兲

Further, to obtain the secular perihelion shift, we recall that Mercury completes 415.19 revolutions per century, so that in conclusion we get, by transforming from radians to seconds,

⌬␻cent =

180 3 ⑀ 2␲ 415.19 3600 = 7.6⬙ , 2n 4n n ␲

共55兲

a value of about 0.3⬙ higher than that obtained with a comprehensive computation,8 while it is found that in 共53兲 the shrinking of the circular component of the orbit is about a couple of hundred meters. In conclusion, we have constructed a general framework to determine the elliptic orbit in function of time that has also the potentiality to be used in some situations in which that elliptical can be considered only a first approximation to the true motion. M. M. D’Eliseo, Can. J. Phys. 85, 1045 共2007兲. M. M. D’Eliseo, Am. J. Phys. 75, 352 共2007兲. T. Needham, Visual Complex Analysis 共Oxford University Press, New York, 1999兲. 4 J. Farlow, J. E. Hall, J. M. McDill, and B. H. West, Differential Equations and Linear Algebra 共Prentice-Hall, Englewood Cliffs, NJ, 2002兲. 5 T. Needham, Visual Complex Analysis 共Oxford University Press, New York, 1999兲, p. 394. 6 L. G. Taff, Celestial Mechanics: A Computational Guide for the Practitioner 共Wiley, New York, 1985兲, pp. 58–66. 7 K. F. Gauss, Determinatio Attractionis Quam in Punctum Quodvis Positonis Datae Exerceret Planeta si Eius Massa per Totam Orbitam Ratione Temporis quo Singulae Partes Describuntur Uniformiter Esset Dispertita 共1818兲 共Werke, Gottingen; Konigliche Gesellschaft der Wissenschaften, 1866兲, Vol. 3, pp. 331–335. 8 M. G. Stewart, Am. J. Phys. 73, 730 共2005兲. 1 2 3

Downloaded 26 Feb 2009 to 192.135.20.186. Redistribution subject to AIP license or copyright; see http://jmp.aip.org/jmp/copyright.jsp