"Problema de los tres cuerpos de Pitágoras": se necesitan algunos puntos de una solución precisa para comparar

nota: si vota a favor (o incluso si no lo hace), no olvide desplazarse hacia abajo y ver también la excelente respuesta: ¡es hermosa!

El problema pitagórico de los tres cuerpos, también conocido como problema de Burrau, es un caso especial del problema general de los tres cuerpos, donde los tres cuerpos tienen masas de 3, 4 y 5, y las condiciones iniciales son tales que comienzan en reposo, en el vértices de un triángulo rectángulo 3-4-5.

He pegado algunas capturas de pantalla de los documentos vinculados aquí .

Puedes ver y leer más en este post

Y mire este video: parece que el tiempo que se muestra en la trama del video es$40\times$tiempo en el papel.

Originalmente, la idea era que podría tener algún significado especial, pero no parece tenerlo. Sin embargo, plantea un gran desafío para los integradores numéricos porque da como resultado varios muy cercanos (~$10^{-4}$) pasa entre pares, y muchos integradores comunes no responderán lo suficientemente rápido con la reducción del tamaño de paso para mantener la precisión numérica.

Esto es lo que me sucedió al usar el integrador ODE predeterminado estándar en SciPy.

Hay algunos trucos para probar dentro de SciPy y, por supuesto, otros integradores disponibles en python, y en realidad puedo implementar algunos métodos de Runge-Kutta de orden superior y escribir mi propio controlador de tamaño de paso hipervigilante . No tiene que ser rápido porque muy pronto, uno de los tres es expulsado y los otros dos se establecen en rotación de dos cuerpos. Esto es bastante común en situaciones de tres cuerpos, en computadoras y en sistemas estelares ternarios que no son lo suficientemente jerárquicos.

Lo que necesito ahora es comparar los resultados con la solución numérica correcta: una tabla con una selección de algunas coordenadas precisas frente al tiempo. Comparar con YouTube no es tan preciso, ¡y tampoco hay garantías de que sean correctos!

¿Alguien sabe dónde puedo encontrar esos números ?

nota: El comentario señala que debo tener cuidado con la palabra "correcto". Estoy buscando resultados usando un solucionador de ODE que funcione bien con ecuaciones rígidas (ver aquí también ) que pueden ser numéricamente inestables , y en este caso se espera que tengan una precisión de, digamos, seis dígitos de precisión por$t=70$.

Aquí hay una salida de muestra y un script. Esto está mal. Puede encontrar buenas soluciones que se muestran en YouTube y otros lugares, pero no puedo encontrar los resultados numéricos para ayudar a mi depuración.

Si desea sugerir una mejora de Python, puede dejar una respuesta o comentario en mi pregunta en 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()