两体系统的Python Euler集成方法没有生成正确的绘图

Question

我一直试图解决一个三体问题，将来应该用于更多的行星，但它不起作用，我应该得到的图是圆形的，但我收到了三体系统的直线。有人知道我如何解决这个问题并得到正确的情节吗？

这是我使用的代码：

import numpy as np
import matplotlib.pyplot as plt
aEarth = 1 #semi-major axis (AU)
eEarth = 0.0167 #eccentricity (no units)
aMercury = 0.387098 #semi-major axis (AU)
eMercury = 0.205635 #eccentricity (no units)

Msun = 1 #Mass of Sun (Solar Mass)
Mearth = 3.0024584*10**(-6) #Mass of Earth (Solar Mass)
Mmercury = 1.65956463*10**(-7) #Mass of Mercury (Solar Mass)
Mes = (Msun + Mearth)
Mms = (Msun + Mmercury)
G = 1 #Gravitational Constant 
apocentreEarth = (aEarth*(1 + eEarth))
apocentreMercury = (aMercury*(1 + eMercury))
vapocentreEarth = (np.sqrt((G*(Mearth+Msun)/aEarth)*((1-eEarth)/(1+eEarth))))
vapocentreMercury = (np.sqrt((G*(Mmercury+Msun)/aMercury)*(1-eMercury/1+eMercury)))
xEarth = apocentreEarth
yEarth = 0.0
zEarth = 0.0
vxEarth = 0.0
vyEarth =(vapocentreEarth)
vzEarth = 0.0
xMercury = apocentreMercury
yMercury = 0.0
zMercury = 0.0
vxMercury = 0.0
vyMercury =(vapocentreMercury)
vzMercury = 0.0
class CelBody(object):
# Constants of nature
def __init__(self, id, name, x0, v0, mass, color, lw):
# Name of the body (string)
self.id = id
self.name = name
# Mass of the body (solar mass)
self.M = mass
# Initial position of the body (au)
self.x0 = np.asarray(x0, dtype=float)
# Position (au). Set to initial value.
self.x = self.x0.copy()
# Initial velocity of the body (au/s)
self.v0 = np.asarray(v0, dtype=float)
# Velocity (au/s). Set to initial value.
self.v = self.v0.copy()
self.a = np.zeros([3], dtype=float)
self.color = color
self.lw = lw
# All Celestial Objects
t = 0
dt = 0.01
Bodies = [
CelBody(0, 'Sun', [0, 0, 0], [0, 0, 0], Msun, 'yellow', 10),
CelBody(1, 'Earth', [xEarth, yEarth, zEarth], [vxEarth, vyEarth, vzEarth], Mearth, 'blue', 3),
CelBody(2, 'Mercury', [xMercury, yMercury, zMercury], [ vxMercury, vyMercury, vzMercury], Mmercury, 'pink', 3),
]
paths = [ [ b.x[:2].copy() ] for b in Bodies]
# loop over ten astronomical years
v = 0
for t in range(0,1000):
# compute forces/accelerations
for body in Bodies:
body.a *= 0
for other in Bodies:
# no force on itself
if (body == other): continue # jump to next loop
rx = body.x - other.x
r3 = (np.sqrt(rx[0]**2+rx[1]**2+rx[2]**2))**3
body.a = -G*other.M*rx/r3
for n, planet in enumerate(Bodies):
# use the Forward Euler algorithm 
planet.a = -G*other.M*rx/r3
planet.v = planet.v + planet.a*dt
planet.x = planet.x + planet.v*dt
paths[n].append( planet.x[:2].copy() )
#print("%10s x:%53s v:%53s"%(planet.name,planet.x, planet.v))

plt.figure(figsize=(8,8))
for n, planet in enumerate(Bodies): 
px, py=np.array(paths[n]).T; 
plt.plot(px, py, color=planet.color, lw=planet.lw)
plt.show()

Answer 1

这些行中存在问题：

planet.v += planet.v + planet.a*dt
planet.x += planet.x + planet.v*dt

它们相当于

planet.v = 2*planet.v + planet.a*dt
planet.x = 2*planet.x + planet.v*dt

这不是你想要的。

要么不使用+=，要么将这些行更改为：

planet.v += planet.a*dt
planet.x += planet.v*dt

仍然存在一个问题：第一行更改planet.v，然后第二行使用新的planet.v来更新planet.x，这不是显式欧拉积分应该如何进行的。(代码abou中有一条注释使用辛欧拉，但对这个答案的注释说你打算使用前向(即显式(欧拉。(

对于这个系统，一个简单的解决方案是切换语句：

planet.x += planet.v*dt
planet.v += planet.a*dt

可能还有其他问题。如果您需要更多帮助，请尽可能简化代码，以创建最小可复制示例。现在，似乎有很多不相关的变量被声明，两个不同的地方是计算欧拉方法的地方，两个地方是分配dt的地方，三个物体是在你说你正在解决两个物体问题时定义的，等等

---

编辑更新后，还剩下几个错误：

• 在开始for body in Bodies:的循环中，计算每个主体的body.a，结果应该累加，因此更新body.a的行应该是

  body.a += -G*other.M*rx/r3
  

  (注意更改为+=(。

• 在应用Euler方法步骤的第二个内部循环(for n, planet in enumerate(Bodies):(中，删除行

  planet.a = -G*other.M*rx/r3
  

  您已经在上一个循环中计算了加速度，并将其存储在planet.a中。

如果我只对现在有问题的代码做了这两个更改，图中显示了您所期望的圆形轨道。