JOURNAL OF

COMPUTATIONAL AND
APPUED MATHEMATICS

ELSEVIER

Journal of Computational and Applied Mathematics 63 (1995) 403-409

Numerical experience with the three-body problem 1

Michal Kfi~.ek *
Mathematical Institute of the Academy of Sciences, ~itnd 25, CZ-11567 Prague 1, Czech Republic
Received 15 September 1994; revised 20 March 1995

Abstract

We discuss the numerical solution of a system of ordinary differential equations that describe the mutual gravitational
influence between three bodies. Several practical examples and some open problems are presented.

Keywords: Runge-Kutta methods; Celestial mechanics; Peculiar orbits

Consider a classical formulation of the three-body problem based on the Newtonian theory, in
which the gravitational forces between bodies depend on their immediate positions and masses, i.e.,
an infinite speed of gravitational influence is assumed. In "nonextremal" cases, such a formulation
yields quite realistic results that are comparable to those obtained from Einstein's theory of relativity
(see [5]).
Let ri = ri(t), i = 1,2, 3, be the so-called radius-vector in 3 describing the position of the ith
body with the mass mi in time t t>0. By Newton's law of gravitation,
Fij = z

mimj(rj -- ri )

Irj

ril3

is the force of the jth body acting on the ith body, where l" I is the Euclidean norm in R 3 and
x = 6.67- 10 -11 kg-lm3s -2 is the gravitational constant. A similar relation holds for the gravitational
force of the ith body on the kth body. Here i , j , k E {1,2,3} and i # j # k ~ i. Moreover, by
Newton's second law, mir~'(t) is the reaction force acting upon the ith body at time t. The threebody problem is thus classically described by the following system of ordinary differential equations
for unknown trajectories rl, r2,r3,
r / ' = x ( m j ( r Z : r-i) m k ( r L : r i ) ~
\ Irj - ri[ 3 q- Irk - ri[ 3 ]

for i = 1,2,3,

(1)

E-mail: krizek@earn.cvut.cz.
t This work was supported by the Grant Agency of the Academy of Sciences of the Czech Republic under grant
No. 119101.
*

0377-0427/95/$09.50 (~) 1995 Elsevier Science B.V. All rights reserved

SSDI 0 3 7 7 - 0 4 2 7 ( 9 5 ) 0 0 0 6 7 - 4

404

M. Kfiek l Journal of Computational and Applied Mathematics 63 (1995) 403-409

where i # j < k # i, and for t = 0 initial positions Pl # P2 # P3

the three bodies are prescribed, i.e.,
ri(O) = p,,

r[(O) = v,

Pt and velocities

for i = 1,2,3.

Vl, [2, V3

of

(2)

Isaac Newton in his Principia (1687) writes: "An exact solution exceeds, if I am not mistaken,
the force of any human mind" (see [12]). Indeed, an analytical solution of (1) and (2) is known
only in a few special cases, in which particular integrals can be computed (see, e.g., [4, 8]). In 1889,
Poincar6 proved (see [9]) that the solution of the three-body problem cannot be explicitly expressed
in terms of known elementary functions, in general. That is why numerical methods are usually
employed.
We observe that the right-hand side of (1) is not Lipschitz continuous (when the ith body is near
to the jth or kth one), which is one of the main sources of numerical troubles. The solution of
(1) and (2) remains uniquely defined (and is regular) as long as all three mutual distances remain
bounded away from zero, i.e., the solution exists locally. However, a global solution of (1) and (2)
need not exist, since the bodies can collide (cf. [6, p. 106]). The aim of this paper is to discuss
numerical solutions of several interesting examples.
Note that stars, planets, asteroids etc. in the three-body problem are replaced by mass points.
We do this because the gravitational field outside a (spherically) homogeneous sphere is equal to
the gravitational field of a mass point that has the same mass and is placed in the centre of the
sphere.

1. Pluto's strange orbit

If [ri(t) - rj(t)[ >> 0 for any i y6 j, then the standard numerical methods for solving initial
problems for ordinary differential equations (see, e.g., [11]) yield quite good results. Such a situation
corresponds, e.g., to the Sun and two planets. In this case we have ml >> mk > 0 for k -- 2, 3. The
choice
Pl = vj = 0,

Pk = (Rk, 0, 0) r, vk = (0, ~ / R k ,

0) 7

for k = 2, 3

(3)

usually produces almost circular orbits. (For the limit case m 2 = m 3 = 0 we would get exactly
circular orbits.) However, if the ratio of periods of the two planets is a fraction of two small
integers, then almost circular orbits may disappear due to resonance. We will demonstrate this in
the following example.
Pluto's orbit is an ellipse (with eccentricity 0.250) which overlaps with Neptune's almost circular
orbit (with eccentricity 0.009). There is an old hypothesis that Pluto was a satellite of Neptune that
escaped in a close encounter with Neptune's largest satellite Triton. Recently another hypothesis
was established. Pluto's orbital period is almost exactly 3 that of Neptune and thus Pluto's elliptical
orbit was caused by this resonance. Note that such a ratio is not valid for any other pair of planets.
According to computations by Malhotra [7], only 5 millions years is enough to achieve Pluto's true
eccentricity of 0.25, starting from the zero eccentricity. This included both analytical and numerical
methods proposed in [13].

M. Kfi~ek / Journal of Computational and Applied Mathematics 63 (1995) 403-409

0.25~ . . . . . . . . . . . . . .

405

true _ _ _

Pluto's.
[
eccentricity ~
.
0.030

computed
.... "---

years

100000
Fig. 1.

We have retested the second hypothesis on the work station HP 720. However, its validity was
not confirmed. We considered the system (1), with the initial conditions (2) and (3), where m~ =
1.99.1030 kg (Sun), m2 = 1.03.1026 kg (Neptune), m 3 ---- 1.4.1022 kg (Pluto), R2 -- 4.4967.1012 m
and by Kepler's third law R 3 = [[
~ j 3 ~2D3
"'2J~1/3 . The gravitational influence of other planets was not
considered. The problem has been solved by the fourth-order Runge-Kutta method in double precision with a fixed time integration step At = 3 months. In this case Pluto's eccentricity was
an oscillating function with a period about 22 500 yr. Its maximum was only 0.035 (see Fig. 1)
and the orbits of Pluto and Neptune did not intersect. Similar results were obtained also for several other fixed integration steps. The perihelion of variable elliptical Pluto's orbit slowly rotated.
The integration over 10 s yr took several days of computer time. Note that the second hypothesis
does not explain why the inclination of Pluto's real orbit to the ecliptic plane is approximately 17
degrees.

2. Motions of comets

When simulating the gravitational influence of a planet with mass m2 on the motion of a comet
with mass m3, we have m~>>m2>>m3,where ml is mass of the Sun. Thus the comet influence upon
the planet and the Sun is almost negligible. A simple version of the corresponding Pascal programme
(which draws trajectories of bodies directly on the screen of usual PC-computers) is available in
[6]. The standard Runge-Kutta methods are very stable even for a fixed time integration step and a
large time interval if [ri(t)- rj(t)[ >>0 for any i ~ j .
The main difficulty in the precise prediction of the trajectory of any comet is not caused by numerical methods but by our inability to establish exactly the initial positions and velocities (see (2)).
Moreover, nongravitational forces affect comet's movement. They come from gas emanation, collisions with meteors and dust, electrical fields etc. These forces are thus unknown. Small changes in
comet velocity then may cause big differences in its trajectory after a long time (cf. [2]). For instance,
in the middle of the year 1993, astronomers calculated that the disintegrated comet Shoemaker-Levy
9 would collide with Jupiter about 23 July, 1994. But we observed collisions of its several pieces
from 16 to 22 July. Therefore, it is almost impossible to predict a clash of a comet with the Earth
many years in advance, even though this is often done. To calculate such a clash more exactly,
mankind needs better observation techniques to establish precisely the initial conditions (2), especially their radial component.

406

M. Kfi~ek l Journal of Computational and Applied Mathematics 63 (1995) 403-409

We give another simple example to demonstrate a high sensitivity of the trajectory with respect
to the initial conditions. Set ml = 1.99. 1030 kg (Sun), m2 = 101 kg (comet) and suppose that the
velocity of a comet in the afelion is v = 5 km s-1 at time t -- 0. Let its distance from the Sun be
p = 1.5.1012 m ( ~ 10 AU) and let T be the associated period. So we have
Pl = vl = 0,

P2 = (P, 0, 0) T,

I)2 = (0, V, 0 ) T,

For t = 0 we shall still consider the perturbed velocity ~ = v + 6v, where 16vl--1 cm s -l only.
Denote by rz the corresponding perturbed trajectory. Then for d ( n ) = I r z ( n T ) - Y2(nT)[, where n
is the number of revolutions, we obtain relatively big distances d ( 1 ) = 1867 kin, d ( 2 ) = 3675 kin,
d(3) = 5625 km etc.

3. Trajectory of a rocket to the Sun

The numerical solution is very sensitive to the initial data (2) and to the choice of time integration step, when ri(t) ~ rj(t) for some t > 0. We often meet this problem in cosmonautics; for
instance, in the computation of a trajectory of a rocket that utilizes Jupiter's gravity to reach the
Sun.
Let us briefly outline why it is convenient to send the rocket to the Sun via Jupiter. The Earth
revolves around the Sun with speed 29 800 m/s. To send the rocket directly to the Sun would mean
giving it the same speed but in the opposite direction. Then the rocket would fall freely directly to
the Sun. Since kinetic energy is proportional to the square of speed, it is much more economical to
send the rocket to Jupiter first. This requires only the second cosmic speed 11 160 m/s. The strong
gravitational field of Jupiter can then direct the rocket towards the Sun, see Fig. 2. The corresponding
demonstration program can be found in [6].
When the rocket is near Jupiter, then their mutual separation is very small compared to the distance
from the Sun. In this case all numerical methods are sensitive to the choice of time integration step,
which should be variable and controlled automatically (see [11]).

start from

the Earth

Fig. 2.

at the
time of start

M. Kfffek l Journal of Computational and Applied Mathematics 63 (1995) 403-409

407

4. Behaviour of a star inside a globular cluster

Globular clusters are objects that contain millions of stars in a relatively small region of about
100 light years in diameter. They exist already 10 ~ yr and it is still not known why they are so
stable (cf. [1]). It has been observed (see, e.g., [3, p. 950]) that the density of stars in a globular
cluster is proportional t o r - a , where r is the distance from cluster's centre of gravity. We assume
that this observation is true for some "reasonable" bounds 0 < ~l < r < 72. Using the standard
spherical coordinates
x = r sin 0 cos q~,
y = r sin 0 sin ~o,
Z ~

r COS 0,

we find that the cluster gravitational potential q/(r) = q/(r, 0, q~) is independent of angles 0 E [0, 7z]
and q> E [0,2rt). Then q / = q/(r) is described by the equation
~32q/

2 0q/

C
in (rl,r2).

- Or
- - i - + r Or -- r4

(4)

The left-hand side corresponds to the Laplace operator in spherical coordinates while the right-hand
side represents the density of gravitational sources (the magnitude of the constant C ~>0 depends
upon the total mass of the cluster). The Eq. (4) can be rewritten as
02

~ r i ( rll ) = -~ .

From here we see that ( ~ / ~ r ) ( r q l ) =

C
C~
q/(r) = ~ r 2 + --r + Co

1 -2
-iCr

+ C 0 and thus
(5)

is a general form of the solution of (4), where Co and C~ are arbitrary constants. The gravitational
force acting on a mass point (star) with mass m is F = - m grad q/. Hence, according to (5) and
Newton's second law, the radius-vector r = r ( t ) of the mass point satisfies (cf. (1))
r

r " = -C7

- C, [r3---~

We again solved this model equation numerically. The trajectory corresponding to C~ > 0 and
C = 0 is an exact ellipse for appropriate initial conditions. For "small" positive C, however, the
computed orbit was nearly elliptic and we observed its perifocus rotation (in the same direction as
that implied by the theory of relativity for C = 0).
To explain the observed rotation of the perifocus, we recall that the gravitational intensity inside
a hollow ball is zero (as the potential is constant there). Hence, the trajectory of a star depends
only on the gravitational forces corresponding to the mass which is inside a ball with radius r (see
Fig. 3). Since r = r ( t ) is variable, the trajectory for C 0 is not an exact ellipse, in general.
To emphasize its shape, we have chosen (in Fig. 3) unrealistic initial conditions, which purposely
produce too large eccentricity.

408

M. KHek/ Journal of Computational and Applied Mathematics 63 (1995) 403-409

Fig. 3.

5. Open problems
Consider now the N-body problem. Its formulation is similar to (1) and (2). We took ml = 100,
m2 . . . . .
m s = 1 for N = 10. The initial conditions of the first N - 1 bodies were of the form (3),
where the Rk were chosen randomly. As all third components are zero, this corresponds to a plane
problem. A nonzero third component was present only in the initial conditions of the Nth body.
This caused a nonzero third component in the other bodies for t > 0, in general. After a very long
time, all the bodies moved chaotically inside a sphere like stars in a globular cluster. The bigger the
mass of the central body, the higher the density of stars with no collisions is possible. Whether the
limit behaviour of the above N bodies tends to a "spherical cluster" is an open problem.
Another nice open problem concerns the structure of Neptune's rings. They mutually intertwine
(see photo taken by Voyager 2 in [10, p. 27]). The mechanism of this peculiar behaviour is not
known. The N-body problem with N ~<20 appeared to be absolutely insufficient to simulate the
intertwining. In the author's opinion, the observed behaviour of the rings cannot be modelled as the
N-body problem for very large N with the advent of supercomputers either. It requires to develop
a completely different model.

Acknowledgements
The author wishes to thank Jana Dafikovfi for valuable suggestions.

References
[1] V.I. Arnold, Small denominators and motion stability problems in classical and celestial mechanics, Uspekhi Mat.
Nauk 18 (1963) 91-192 (in Russian).
[2] Z. Ceplecha, Influx of interplanetary bodies onto Earth, Astronom. Astrophys. 263 (1992) 361-366.
[3] V. Guth, F. Link, J. M. Mohr and B. Sternberk, Astronomie H (Academia, Prague, 1954).
[41 W.W. Heinrich, On new particular integrals of the satellite problem of three bodies, Vestnik Krdl. Ces. Spol. Nauk
IV (1952) 1-45.

M. Kfi~ek l Journal of Computational and Applied Mathematics 63 (1995) 403-409

409

V. Hlavat~, Geometry of Einstein's Unified Field Theory (Noordhoff, Holland, 1958).

M. IZd-i~ek,On the three-body problem, Rozhledy mat.-fyz. 70 (1992) 105-112 (in Czech).
R. Malhotra, The origin of Pluto's peculiar orbit, Nature 365 (1993) 819-821.
C. Marchal, The Three-Body Problem (Elsevier, Amsterdam, 1991).
H. Poincar~, Sur le probl~me des trois corps et les ~quations de ia dynamique Acta Math. 13 (1890) 1-270.
B. Sicardy, Les anneaux de Neptune, La Recherche, nr. 261, 25 (1994) 22-29.
E. Vitfisek, Principles of the Theory of Numerical Methods for Solving Differential Equations (Academia, Prague,
1994).
[12] C. Wilson, The three-body problem, in: I. Grattan-Guinness, Ed., Companion Encyclopedia of the History and
Philosophy of the Mathematical Sciences, 2 (Routledge, London, New York, 1994) 1054-1062.
[13] J. Wisdom and M. Holman, Symplectic maps for the N-body problem, Astronom. J. 102 (1991) 1528-1538.

[5]
[6]
[7]
[8]
[9]
[10]
[11]