следующий шаг в тестировании решателя ОДУ для пифагорейской задачи трех тел

Я столкнулся с трудностями, пытаясь интегрировать Задачу Пифагора о трех телах с помощью scipy.odeint. После небольшого осмотра и поиска в Интернете я нашел следующее в этой очень интересной интеграции обсуждение/руководство:

"После обсуждения масштабирования единиц измерения в следующем разделе в следующих разделах описано множество различных алгоритмов интегрирования. Автор рекомендует после написания собственной программы интегрирования в соответствии с одним из этих алгоритмов начать упражнения по интегрированию с рисунка «восьмерка», так как ее легко интегрировать из-за ее стабильности и того факта, что близких столкновений вообще не происходит. После этого вы можете попытаться решить задачу Пифагора. Задача Пифагора трудно интегрируется. A необходимо использовать очень точный интегратор, способный справиться с многочисленными близкими столкновениями."

Итак, мой основной вопрос: есть ли другие библиотеки ODE для Python, на которые я мог бы обратить внимание, в соответствии с приведенным выше предложением? В качестве альтернативы, может ли кто-нибудь помочь мне понять, как уговорить odeint работать здесь? scipy.odeint всегда «просто работал» прямо из коробки всякий раз, когда я его использовал, поэтому на этот раз я был удивлен.

В этом видео и в этом видео

примечание: заголовок не является опечаткой - бот блокирует слово "проблема" в заголовке.

Я собираюсь опубликовать свою первую попытку реализации ниже. Буду рад комментариям, как лучше написать. Регулируя tol (а иногда и интервал в t, что странно, потому что это интерполяция, а не фактические временные шаги для scipy.odeint). Как только я смог создать правильно выглядящий сюжет (вы можете увидеть их повсюду на интернет), но не помню как.

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