note: if you up-vote (or even if you don't), don't forget to scroll down and see the excellent answer as well - it's beautiful!

The Pythagorean Three Body Problem also know as Burrau's problem is a special case of the general three body problem, where the the three bodies have masses of 3, 4, and 5, and the initial conditions are such that they begin at rest, at the vertices of a 3-4-5 right triangle.

I've pasted some screenshots from the papers linked here.

You can see and read more in this posting

And watch this video - it looks like time displayed in the plot in the video is $40\times$ time in the paper.

The idea was originally that it might hold some special significance, but it doesn't seem to. However, it does pose a big challenge to numerical integrators because it results in several very close (~$10^{-4}$) passes between pairs, and many common integrators will not respond quick enough with step-size reduction to maintain numerical accuracy.

This is what has happened to me using the standard default ODE integrator in SciPy.

There are some tricks to try within SciPy, and of course other integrators available in python, and actually I can just implement some higher order Runge-Kutta methods and write my own hyper-vigilant step size handler. It does't have to be fast because fairly soon, one of the three is ejected and the other two settle down to two-body rotation. This is pretty common in three body situations, in computers and in ternary star systems that aren't sufficiently hierarchical.

What I do need now is to compare results to the correct numerical solution - a table with a selection of some precise coordinates vs time. Comparing to YouTube isn't so accurate, and no guarantees those are right either!

Does anyone know where I can find such numbers?

note: The comment points out I should be careful with the word "correct." I'm looking for results using an ODE solver that works well with stiff equations (see here also) that may be numerically unstable, and in this case are expected to be accurate to - say - six digits of accuracy by $t=70$.

Here is a sample output and script. This is wrong. You can find nice solutions displayed in YouTube and other places, but I can't find the numerical results to help my debugging.

If you want to suggest python improvement, you can leave an answer or comment at my question in stackoverflow

Wrong Answer

def deriv(X, t):

    Y[:6] = X[6:]

    r34, r35, r45 = X[2:4]-X[0:2], X[4:6]-X[0:2], X[4:6]-X[2:4]
    thing34 = ((r34**2).sum())**-1.5
    thing35 = ((r35**2).sum())**-1.5
    thing45 = ((r45**2).sum())**-1.5

    Y[6:8]   =  r34*thing34*m4 + r35*thing35*m5
    Y[8:10]  =  r45*thing45*m5 - r34*thing34*m3
    Y[10:12] = -r35*thing35*m3 - r45*thing45*m4

    return Y

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint as ODEint

# Pythagorean Three Body Problem
# This script WILL NOT solve it yet, just for illustration of the problem

m3, m4, m5 = 3.0, 4.0, 5.0

x0 = [1.0, 3.0] + [-2.0, -1.0] + [1.0, -1.0]
v0 = [0.0, 0.0] + [ 0.0,  0.0] + [0.0,  0.0] 
X0 = np.array(x0 + v0)

t = np.linspace(0, 60,  50001)

Y = np.zeros_like(X0)

tol  = 1E-9 # with default method higher precision causes failure
hmax = 1E-04
answer, info = ODEint(deriv, X0, t, rtol=tol, atol=tol,
                      hmax=hmax, full_output=True)

xy3, xy4, xy5 = answer.T[:6].reshape(3,2,-1)
paths         = [xy3, xy4, xy5]

plt.subplot(2, 1, 1)
for x, y in paths:
    plt.plot(x, y)
for x, y in paths:
    plt.plot(x[:1], y[:1], 'ok')
plt.xlim(-6, 6)
plt.ylim(-4, 4)
plt.title("This result is WRONG!", fontsize=16)
for x, y in paths:
    plt.plot(t, x)
plt.ylim(-6, 4)
for x, y in paths:
    plt.plot(t, y)
plt.ylim(-6, 4)

enter image description here enter image description here enter image description here enter image description here enter image description here enter image description here enter image description here enter image description here

  • 1
    $\begingroup$ Unless I'm completely wrong, any 3-body simulation with close approaches is going to be chaotic, that is to say sensitively dependent on initial conditions. You'll be at the mercy of not only time quantization but of floating point numerical artifacts; there is no "correct answer" in a finite/discrete system to compare to. $\endgroup$ May 21, 2016 at 18:32
  • $\begingroup$ @RussellBorogove Thanks, I'll add something about the "c" word (correct). it's true that chaotic systems are very sensitive to initial conditions, and numerical artifacts limit accuracy, but those don't mean the problem is stochastic or non-deterministic. Two programmers with two different libraries on two different computers using double precision (8 bytes) should be able to solve this problem certainly to - say - six digits of accuracy - to the point ($t>60$) where the $m=3$ object is ejected. There's no uncertainty principle or time quantization here, just solving a series of 4 ODEs. $\endgroup$
    – uhoh
    May 22, 2016 at 0:04
  • $\begingroup$ Floating point operations on standard CPUs can yield different results for operations as similar as adding three numbers in different orders. Rounding effects can give different results for (3 steps of delta-t = 2) versus (2 steps of delta-t = 3). As soon as you get a divergence, the divergence will start to amplify. $\endgroup$ May 22, 2016 at 0:29
  • $\begingroup$ @ that's right. I'm starting with 16 digit numbers, hoping for about 6 digits of accuracy at the end ($t=70$). Right now I have tol set large $10^{-10}$ because the default routine isn't able to handle the problem. See this question for more examples - problem D5 in particular. Of course different computers can give answers that can be slightly different, and round-off errors can lead to much bigger errors with this kind of problem (chaotic). $\endgroup$
    – uhoh
    May 22, 2016 at 0:58
  • $\begingroup$ @RussellBorogove I'm looking for a solution where someone has used "numerical overkill" (lots of steps, stiff-friendly method, maybe quad precision) to get a good answer, to which I can compare. $\endgroup$
    – uhoh
    May 22, 2016 at 0:58

1 Answer 1


I just ran it, and mine look pretty much like those in the paper.

See some coordinates at the bottom.

0-10 10-20 20-30 30-40 40-50 50-60 60-70

Here are some {x,y} coordinates at the times in the left column:

0.      {1.,3.}                 {-2.,-1.}               {1.,-1.}
5.      {2.46917,-1.22782}      {-2.2782,-0.20545}      {0.34106,0.901049}
10.     {0.77848,0.141392}      {-2.02509,0.0972194}    {1.15299,-0.162611}
15.     {1.41845,0.686214}      {-2.00654,0.0599408}    {0.754159,-0.459681}
20.     {3.00429,0.511925}      {-1.38863,-0.470476}    {-0.691674,0.0692257}
25.     {2.2699,-0.0832}        {-2.63692,-0.426417}    {0.747596,0.391054}
30.     {0.85634,2.28709}       {-0.877984,-0.865964}   {0.188583,-0.679485}
35.     {0.0273748,0.895529}    {0.942553,-1.60223}     {-0.770468,0.744467}
40.     {-0.622004,1.85832}     {0.173545,-2.36841}     {0.234367,0.779737}
45.     {-0.657058,2.53557}     {1.61355,-1.23947}      {-0.896608,-0.529771}
50.     {-2.70146,-3.79723}     {1.50595,0.960811}      {0.416122,1.50969}
55.     {-2.75171,-4.29907}     {1.72673,0.97731}       {0.269648,1.7976}
60.     {0.743681,1.93961}      {0.263967,-0.731477}    {-0.657382,-0.578586}
65.     {4.05348,11.7131}       {-1.0722,-3.92197}      {-1.57432,-3.8903}
70.     {6.93108,20.2566}       {-1.99418,-6.87252}     {-2.5633,-6.65594}

That was all with 30 digits working precision. Checking the final total energy and total angular momentum against the initial conditions, with 30 working digits the results are good to 10 digits. With 50 working digits, the results are good to 20 digits. With machine precision (about 15 working digits), the results are good to five to six digits, which is still pretty good considering the close approaches.

  • $\begingroup$ Wow - beautiful! I had a feeling you'd come through :) How about a 17x6 table of numbers - 3 pairs of (x, y) for t = 0, 5, 10... 70 - uneven time spacing is fine. Also, maybe a quick check that your solution seems stable against changes in whichever tolerance parameters you might have. It's not "proof of accuracy", but it's good enough for me. Thanks very much!! $\endgroup$
    – uhoh
    May 22, 2016 at 2:22
  • 1
    $\begingroup$ I ran it with a higher working precision (30 decimal digits) and got the same results. $\endgroup$
    – Mark Adler
    May 22, 2016 at 2:43
  • 1
    $\begingroup$ Mathematica, right? I'm curious, did you use the default settings or some custom ones? $\endgroup$ May 23, 2016 at 13:53
  • 1
    $\begingroup$ Yes. The plot colors should be a dead giveaway. I used NDSolve with InterpolationOrder -> All, WorkingPrecision -> 30, MaxSteps -> 10^5. $\endgroup$
    – Mark Adler
    May 23, 2016 at 15:26
  • 2
    $\begingroup$ There is a SymplecticPartitionedRungeKutta option available, but I didn't use it. I used the default methods, which choose a predictor-corrector and a backward differentiation method, depending on stiffness. Then the final total energy really is a good measure of the result quality, since there is nothing explicit in the integration method, other than the equations of motion, that would assure its conservation. $\endgroup$
    – Mark Adler
    May 23, 2016 at 16:50

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.