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
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.figure()
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)
plt.subplot(4,1,3)
for x, y in paths:
plt.plot(t, x)
plt.ylim(-6, 4)
plt.subplot(4,1,4)
for x, y in paths:
plt.plot(t, y)
plt.ylim(-6, 4)
plt.show()
tolset 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$