<p>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 <a href="https://arxiv.org/ftp/arxiv/papers/1605/1605.05317.pdf" rel="noreferrer">Existence and Uniqueness theorems of differential equations</a>. You'd be better to ask questions like this on the Mathematics stack exchange.</p>

<p>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.</p>

<p>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 <a href="https://astronomy.stackexchange.com/questions/23239/is-the-three-body-system-unique#comment39913_23239">userLTK states</a>, you can write down approximate solutions in the <a href="http://scienceworld.wolfram.com/physics/RestrictedThree-BodyProblem.html" rel="noreferrer">restricted 3-body problem</a>, where one mass is significantly less than the other two and orbits under specific conditions.</p>

<p>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.</p>

<p>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.</p>

<pre><code>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),
               1E11*np.sin(210*np.pi/180)])
v2 = np.array([3E4*np.cos(300*np.pi/180),
               3E4*np.sin(300*np.pi/180)])

# Define body 3 parameters
m3 = 2E30
x3 = np.array([1E11*np.cos(330*np.pi/180),
               1E11*np.sin(330*np.pi/180)])
v3 = np.array([3E4*np.cos(60*np.pi/180),
               3E4*np.sin(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):
        a.append(0)
        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
        vHist[i].append(list(vi))

    # Update positions
    for i,xi,vi in zip(range(n),x,v):
        xi += vi * dt
        xHist[i].append(list(xi))

    # Update time and end simulation if past tEnd
    t += dt
    if dt &gt; 0 and t &gt; tEnd: break
    if dt &lt; 0 and t &lt; 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
figure()
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)
gca().set_aspect(1)
gca().set_xticks([])
gca().set_yticks([])
show(block = False)
</code></pre>

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

<p><strong>Path of three masses for $t&lt;t_0$</strong>
<a href="https://i.sstatic.net/RvlD3.png" rel="noreferrer"><img src="https://i.sstatic.net/RvlD3.png" alt="enter image description here"></a>
In this scenario, the three masses end up in the contrived scenario. I ran the simulation backwards by setting the timestep <code>dt</code> and the end time <code>tEnd</code> to be negative.</p>

<p><strong>Path of three masses for $t_0&lt;t$</strong>
<a href="https://i.sstatic.net/Tc5Mc.png" rel="noreferrer"><img src="https://i.sstatic.net/Tc5Mc.png" alt="enter image description here"></a>
From here, the simulation is run forwards with positive <code>dt</code> and <code>tEnd</code>, starting from the same contrived scenario as above.</p>

<p>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.</p>

<p><strong>Contrived Scenario with $m_1=m_2=m_3$, $|v_1|=|v_2|=|v_3|$, and all positions $120^o$ from each other</strong>
<a href="https://i.sstatic.net/fYl6M.png" rel="noreferrer"><img src="https://i.sstatic.net/fYl6M.png" alt="enter image description here"></a>
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 <em>any</em> 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.</p>

CC BY-SA 3.0
2017-10-21T06:12:19.870Z
Thank you for your comprehensive answer! Just out of curiosity, how long did you run your experiment until the first system ejects one body?
xtt
CC BY-SA 3.0
2017-10-21T13:09:48.707Z
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.
zephyr
CC BY-SA 3.0
2017-10-23T13:57:50.410Z
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](https://en.wikipedia.org/wiki/Leapfrog_integration); another popular choice of symplectic integrator is Verlet integration, as mentioned on that Wikipedia page.
PM 2Ring
CC BY-SA 3.0
2017-10-23T14:01:01.823Z
@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.
zephyr
CC BY-SA 3.0
2017-10-20T15:08:27.190Z
2017-10-20T15:08:27.190Z
zephyr
<p>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}&gt;t$. But is it possible to determine a unique state in time $\bar{t}$ where $\bar{t}&lt;t$? </p>

<p>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?</p>

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

CC BY-SA 3.0
2017-10-20T11:19:43.190Z
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.
StephenG - Help Ukraine
CC BY-SA 3.0
2017-10-20T11:56:44.533Z
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.
userLTK
CC BY-SA 3.0
2017-10-20T13:42:19.707Z
Ask Cixin Liu :-)
Carl Witthoft
<orbit><stellar-dynamics>
Is the three-body system "unique"?
xtt
# Is the three-body system "unique"?

(Asked by: username_0 on 2017-10-20)

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?

**Tags:** orbit, stellar-dynamics

> **Comments:**
> **username_1**: 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.  
> **username_2**: 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.  
> **username_3**: Ask Cixin Liu :-)  

---

## 1 answers

### Answer by username_4 on 2017-10-20. Score: 8 | Accepted answer

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](https://arxiv.org/ftp/arxiv/papers/1605/1605.05317.pdf). 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 [username_2 states](https://astronomy.stackexchange.com/questions/23239/is-the-three-body-system-unique#comment39913_23239), you can write down approximate solutions in the [restricted 3-body problem](http://scienceworld.wolfram.com/physics/RestrictedThree-BodyProblem.html), 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),
               1E11*np.sin(210*np.pi/180)])
v2 = np.array([3E4*np.cos(300*np.pi/180),
               3E4*np.sin(300*np.pi/180)])

# Define body 3 parameters
m3 = 2E30
x3 = np.array([1E11*np.cos(330*np.pi/180),
               1E11*np.sin(330*np.pi/180)])
v3 = np.array([3E4*np.cos(60*np.pi/180),
               3E4*np.sin(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):
        a.append(0)
        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
        vHist[i].append(list(vi))

    # Update positions
    for i,xi,vi in zip(range(n),x,v):
        xi += vi * dt
        xHist[i].append(list(xi))

    # 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
figure()
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)
gca().set_aspect(1)
gca().set_xticks([])
gca().set_yticks([])
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](https://i.sstatic.net/RvlD3.png)](https://i.sstatic.net/RvlD3.png)
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](https://i.sstatic.net/Tc5Mc.png)](https://i.sstatic.net/Tc5Mc.png)
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](https://i.sstatic.net/fYl6M.png)](https://i.sstatic.net/fYl6M.png)
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.

> **Comments:**
> **username_0**: Thank you for your comprehensive answer! Just out of curiosity, how long did you run your experiment until the first system ejects one body?  
> **username_4**: 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.  
> **username_5**: 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](https://en.wikipedia.org/wiki/Leapfrog_integration); another popular choice of symplectic integrator is Verlet integration, as mentioned on that Wikipedia page.  
> **username_4**: @username_5 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.  

---

