"Pythagoräisches Dreikörperproblem" - brauche einige Punkte aus einer genauen Lösung zum Vergleich

Hinweis: Wenn Sie positiv abstimmen (oder auch nicht), vergessen Sie nicht, nach unten zu scrollen und sich auch die hervorragende Antwort anzusehen - es ist wunderschön!

Das pythagoräische Drei-Körper-Problem , auch bekannt als Burraus Problem, ist ein Sonderfall des allgemeinen Drei-Körper-Problems, bei dem die drei Körper Massen von 3, 4 und 5 haben und die Anfangsbedingungen so sind, dass sie im Ruhezustand beginnen Eckpunkte eines rechtwinkligen 3-4-5-Dreiecks.

Ich habe einige Screenshots aus den hier verlinkten Papieren eingefügt .

In diesem Beitrag können Sie mehr sehen und lesen

Und sehen Sie sich dieses Video an - es sieht so aus, als ob die Zeit in der Handlung im Video angezeigt wird$40\times$Zeit in der Zeitung.

Die Idee war ursprünglich, dass es eine besondere Bedeutung haben könnte, aber es scheint nicht so zu sein. Es stellt jedoch eine große Herausforderung für numerische Integratoren dar, da es zu mehreren sehr nahen (~$10^{-4}$) zwischen Paaren durchläuft, und viele gängige Integratoren reagieren nicht schnell genug mit einer Reduzierung der Schrittgröße, um die numerische Genauigkeit aufrechtzuerhalten.

Dies ist mir mit dem Standard-ODE-Integrator in SciPy passiert.

Es gibt einige Tricks, die Sie in SciPy ausprobieren können, und natürlich andere Integratoren, die in Python verfügbar sind, und tatsächlich kann ich einfach einige Runge-Kutta-Methoden höherer Ordnung implementieren und meinen eigenen hyper-wachsamen Step-Size-Handler schreiben . Es muss nicht schnell sein, denn ziemlich bald wird einer der drei ausgeworfen und die anderen beiden beruhigen sich auf eine Zwei-Körper-Rotation. Dies ist ziemlich üblich in Drei-Körper-Situationen, in Computern und in ternären Sternensystemen, die nicht ausreichend hierarchisch sind.

Was ich jetzt brauche, ist, die Ergebnisse mit der richtigen numerischen Lösung zu vergleichen - eine Tabelle mit einer Auswahl einiger genauer Koordinaten gegenüber der Zeit. Der Vergleich mit YouTube ist nicht so genau, und es gibt auch keine Garantie dafür, dass diese richtig sind!

Weiß jemand wo ich solche Nummern finden kann ?

Hinweis: Der Kommentar weist darauf hin, dass ich mit dem Wort "richtig" vorsichtig sein sollte. Ich suche nach Ergebnissen mit einem ODE-Löser, der gut mit steifen Gleichungen funktioniert (siehe auch hier ), die möglicherweise numerisch instabil sind und in diesem Fall voraussichtlich eine Genauigkeit von - sagen wir - sechs Stellen haben$t=70$.

Hier ist eine Beispielausgabe und ein Skript. Das ist falsch. Sie können nette Lösungen finden, die in YouTube und anderen Orten angezeigt werden, aber ich kann die numerischen Ergebnisse nicht finden, um mein Debugging zu unterstützen.

Wenn Sie Python-Verbesserungen vorschlagen möchten, können Sie eine Antwort oder einen Kommentar zu meiner Frage in Stackoverflow hinterlassen

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()