Given a state of an ideal 3 body system (i.e., without external interference) in time $t$: the velocity $v_{i,t}$, mass $m_{i,t}$ and position $x_{i,t}$ for $i\in \{1,2,3\}$, using numerical method it is possible to determine any state in time $\hat{t}$, where $\hat{t}>t$. But is it possible to determine a unique state in time $\bar{t}$ where $\bar{t}<t$?

In another word, is it sufficient to know a single state of the system for any given time to deduce the states of the system for all time?

Equivalently, will two different states of the 3-body system results in the same state (maybe at different time) in the future?

  • 1
    $\begingroup$ You can never have perfect knowledge of the state of a system, so there would always be uncertainty in the initial state. Also the 3 body problem has no general analytic solution, so we would have to approximate with numerical methods, which essentially guarantees additional errors. $\endgroup$ Oct 20, 2017 at 11:19
  • 1
    $\begingroup$ There are, also, two kinds of 3 body problems. Long term stable and Chaotic. The Sun-Earth-Moon or Sun-Earth-Mars are examples of long term stable 3 body problems where the orbits, but not necessarily the exact positions, are long term predictable. The chaotic 3-body problem is usually not long term stable. Both can be approximated quite far in the future with computers and some tedious mathematical programing, but as Stephen G points out, errors multiply the further out you go. $\endgroup$
    – userLTK
    Oct 20, 2017 at 11:56
  • 1
    $\begingroup$ Ask Cixin Liu :-) $\endgroup$ Oct 20, 2017 at 13:42

1 Answer 1


What you're really asking about is less to do with Astronomy and more to do with mathematics. You're basically asking if, given a system of differential equations, will a unique solution exist for all time? For an answer, you should check out the Existence and Uniqueness theorems of differential equations. You'd be better to ask questions like this on the Mathematics stack exchange.

However, to discuss the particular, astronomical case you've asked about, the answer is yes, you can run that system both forwards and backwards if you know some initial state. Newtonian mechanics is completely deterministic in that if you know all the equations of motion involved, as well as the entire state of the system at a given time, you can figure out the state of that system at any other time, both in the past and the future.

To speak to your particular 3-body orbiting problem though, I'll say that the system of equations is not solvable in closed form - that is, you can't write down an analytic solution to the equations like you could for the 2-body case. As userLTK states, you can write down approximate solutions in the restricted 3-body problem, where one mass is significantly less than the other two and orbits under specific conditions.

To get a solution at any time $\hat{t}$, you need to use numerical methods. Of course numerical methods are inherently flawed. Numerical errors build up the longer you simulate due to time steps which are not infinitesimal, from errors within the numerical algorithms used, and general floating point errors. In theory, if you had a computer with infinite precision and infinite computing power, you could solve a 3-body (or n-body) system perfectly, but we live in the real world where such things are impossible.

To prove though, that you can figure out the state at any time in the past or present for a three body body, I've written a basic simulation in Python 3. It can run both forwards and backwards from a given start condition and start time. Essentially it puts three nearly identical masses in contrived starting positions and velocities. Below the code are plots of results.

import numpy as np
from numpy.linalg import norm
from matplotlib.pyplot import *
from time import time

# Define physical constants
G = 6.67408E-11          # Gravitational Constant, m^3 kg^-1 s^-2

# Define body 1 parameters
m1 = 2.2E30                # Mass, kg
x1 = np.array([0,1E11])  # Position, m
v1 = np.array([-3.5E4,0])  # Velocity, m/s

# Define body 2 parameters
m2 = 1.9E30
x2 = np.array([1E11*np.cos(210*np.pi/180),
v2 = np.array([3E4*np.cos(300*np.pi/180),

# Define body 3 parameters
m3 = 2E30
x3 = np.array([1E11*np.cos(330*np.pi/180),
v3 = np.array([3E4*np.cos(60*np.pi/180),

# Define simulation parameters
n = 3                    # Number of bodies, unitless
t = 0                    # Simulation time, s
dt = 1E4                 # Simulation time step, s
tEnd = 1E8               # Simulation end time, s
m = np.array((m1,m2,m3)) # All masses
x = np.vstack((x1,x2,x3))# All positions
v = np.vstack((v1,v2,v3))# All velocities
xHist = [[list(x1)],[list(x2)],[list(x3)]]
vHist = [[list(v1)],[list(v2)],[list(v3)]]

# Simulate until end time is reached
start = time()
while True:
    # Calculate acceleration
    a = []
    for i in range(n):
        for j in range(n):
            if i == j: continue
            a[-1] += - G * m[j] / norm(x[i]-x[j])**3 * (x[i] - x[j])

    # Update velocities
    for i,vi,ai in zip(range(n),v,a):
        vi += ai * dt

    # Update positions
    for i,xi,vi in zip(range(n),x,v):
        xi += vi * dt

    # Update time and end simulation if past tEnd
    t += dt
    if dt > 0 and t > tEnd: break
    if dt < 0 and t < tEnd: break

end = time()
print('Simulation finished in {:.4f} seconds.'.format(end-start))

# Convert xHist and vHist to np arrays
xHist = np.array(xHist)
vHist = np.array(vHist)

# Plot everything up
for i,c in enumerate(['or','ob','oc']):
    # Plot starting positions
    plot(xHist[i,0,0], xHist[i,0,1], c)
for i,c in enumerate(['-r','-b','-c']):
    # Plot path of star
    plot(xHist[i,:,0], xHist[i,:,1], c)
show(block = False)

Note, the plots show the initial positions of the stars as the points and then trace out their paths over time.

Path of three masses for $t<t_0$ enter image description here In this scenario, the three masses end up in the contrived scenario. I ran the simulation backwards by setting the timestep dt and the end time tEnd to be negative.

Path of three masses for $t_0<t$ enter image description here From here, the simulation is run forwards with positive dt and tEnd, starting from the same contrived scenario as above.

Note how chaotic and unstable this system is. The entire system only remains a 3-body system for less than 6 years. Before that the three masses are separate and doing their own thing. They "coincidentally" meet (because I set it up so they should), orbit around each other for a little less than 6 years, and one gets ejected, resulting in the other two continuing to orbit one another.

Contrived Scenario with $m_1=m_2=m_3$, $|v_1|=|v_2|=|v_3|$, and all positions $120^o$ from each other enter image description here Just for fun, a "stable" 3-body problem with all stars of equal mass and orbits. This is really an unstable equilibrium orbit and any perturbations will screw it up and you'll see what you saw in the above two graphs. In fact, if you run this long enough, the numerical instabilities of my code will result in the orbits breaking down. This numerical instability, as I said above, is inherent in any numerical solution and cannot be overcome, only minimized. I find that, using the numerical method I did, my system is resistant to numerical instabilities for about 7 years. If I want to run this any longer (and continue to be accurate), I need more robust numerical methods.

  • 1
    $\begingroup$ Thank you for your comprehensive answer! Just out of curiosity, how long did you run your experiment until the first system ejects one body? $\endgroup$
    – xtt
    Oct 21, 2017 at 6:12
  • 1
    $\begingroup$ As I said above, the entire time was about 6 years. It was about evenly spilt, with the $t<t_0$ taking about three years and the $t_0<t$ taking about three years. Overall, not a very stable system. $\endgroup$
    – zephyr
    Oct 21, 2017 at 13:09
  • 2
    $\begingroup$ FWIW, you can reduce the numerical errors by using a better integration technique than Euler integration (which is what your code uses). For orbital simulations you should use a symplectic integrator (i.e., one that conserves energy). My preference is for the synchronized form of Leapfrog integration; another popular choice of symplectic integrator is Verlet integration, as mentioned on that Wikipedia page. $\endgroup$
    – PM 2Ring
    Oct 23, 2017 at 13:57
  • 2
    $\begingroup$ @PM2Ring That's definitely true. I wasn't trying to make a very robust simulation. Just something quick and dirty to prove a point. But yes, a real simulation would be better served by using something other than the basic Euler integration. $\endgroup$
    – zephyr
    Oct 23, 2017 at 14:01

You must log in to answer this question.

Not the answer you're looking for? Browse other questions tagged .