Python 两个 body 系统的欧拉积分方法未生成正确的绘图

Question

我一直在尝试解决两个 body 问题，将来应该用于更多行星，但它没有用，我应该得到的情节是循环的，但我我收到了两个 body 系统的直线。有谁知道我该如何解决这个问题并获得正确的情节？

这是我使用的代码：

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，然后第二行用newplanet.v更新planet.x ，这不是显式欧拉积分应该如何进行。（代码中有一条关于使用辛欧拉的评论，但是对此答案的评论说你打算使用前向（即显式）欧拉。）

对于这个系统，一个简单的修复是切换语句：

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

可能还有其他问题。如果您需要更多帮助，请尽可能简化代码以创建minimal reproducible example。就像现在一样，看起来有很多不相关的变量声明，两个不同的地方你计算欧拉的方法，两个地方你分配 dt，当你说你正在解决两个 body问题等

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

在开始 for body in Bodies: 的循环中，您在其中为每个 body 计算 body.a，结果应该累加，因此更新 body.a 的行=22=] 应该是
```
            body.a += -G*other.M*rx/r3
```
（注意 += 的变化）。
在第二个内部循环 (for n, planet in enumerate(Bodies):) 中，您在其中应用欧拉方法步骤，删除行
```
        planet.a = -G*other.M*rx/r3
```
您已经计算了上一个循环中的加速度并将它们存储在 planet.a 中。

如果我只对现在问题中的代码进行这两项更改，绘图会显示您预期的圆形轨道。

Python 两个 body 系统的欧拉积分方法未生成正确的绘图

Python Euler Integration Method for two body system not producing correct plot

python

integration

numpy

overflow