matlab ode45三体问题,“毕达哥拉斯3Body Proxblem”ODE解算器测试的下一步

我在尝试使用scipy.odeint整合毕达哥拉斯三体问题时遇到了困难。经过一番检查和网络搜索，我在这个非常有趣的集成中发现了以下内容discussion/tutorial："After a discussion of the unit scaling in the next section many different integration algorithms are described in following sections. The author recommends, after writing your own integration program according to one of these algorithms, to begin the integration exercises with the figure “eight”, since it is easy to integrate due to its stability and the fact, that close encounters do not occur at all. Afterwards, you may try to solve the Pythagorean problem. The Pythagorean problem is difficult to integrate. A very accurate integrator must be used which is able to cope with the numerous close encounters."

因此，我的主要问题是：是否有其他python-ODE库可以参照上面的建议？或者，有人能帮我理解如何诱使odeint在这里工作吗？scipy.odeint每当我使用它时，它总是“刚刚工作”，所以这次我很惊讶。在

这个video和这个video都有很好的模拟结果

注意：标题不是打字错误-标题中有一个bot阻止了单词“problem”。在

我将在下面发布我的第一次尝试实现。我欢迎评论如何写得更好。通过调整tol(有时还有t中的间距，这很奇怪，因为这是插值，而不是{}的实际时间步长)。有一次我能画出一个看起来很正确的图(你可以在internet上看到它们)，但我不记得是怎么画的。在

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