两体问题中的 Python Euler 方法实现不起作用

Question

我目前正在尝试解决二体问题，然后我可以升级到更多行星，但它不起作用。它正在输出我不可能的位置。有谁知道是什么原因造成的？

这是我使用的代码:

day = 60*60*24
# Constants
G = 6.67408e-11
dt = 0.1*day
au = 1.496e11
t = 0


class CelBody:

    def __init__(self, id, name, x0, y0, z0, vx0, vy0, vz0, mass, vector, ax0, ay0, az0, totalforcex, totalforcey, totalforcez):
        self.ax0 = ax0
        self.ay0 = ay0
        self.az0 = az0

        self.ax = self.ax0
        self.ay = self.ay0
        self.az = self.az0

        # Constants of nature
        # Universal constant of gravitation
        self.G = 6.67408e-11
        # Name of the body (string)
        self.id = id
        self.name = name
        # Initial position of the body (au)
        self.x0 = x0
        self.y0 = y0
        self.z0 = z0
        # Position (au). Set to initial value.
        self.x = self.x0
        self.y = self.y0
        self.z = self.z0
        # Initial velocity of the body (au/s)
        self.vx0 = vx0
        self.vy0 = vy0
        self.vz0 = vz0
        # Velocity (au/s). Set to initial value.
        self.vx = self.vx0
        self.vy = self.vy0
        self.vz = self.vz0
        # Mass of the body (kg)
        self.M = mass
        # Short name
        self.vector = vector

        self.totalforcex = totalforcex
        self.totalforcey = totalforcey
        self.totalforcez = totalforcez

# All Celestial Bodies

forcex = 0
forcey = 0
forcez = 0

Bodies = [
    CelBody(0, 'Sun', 1, 1, 1, 0, 0, 0, 1.989e30, 'sun', 0, 0, 0, 0, 0, 0),
    CelBody(1, 'Mercury', 1*au, 1, 1, 0, 29780, 0, 3.3e23, 'earth', 0, 0, 0, 0, 0, 0),
    ]

leftover_bin = []
templistx = []
templisty = []
templistz = []

for v in range(365242):
    for n in range(len(Bodies)):
        #Need to initialize the bodies

        planetinit = Bodies[n]

        for x in range(len(Bodies)):
            # Temporary lists and initial conditions
            planet = Bodies[x]

            if (planet == planetinit):
                pass

            else:
                rx = Bodies[x].x - Bodies[n].x
                ry = Bodies[x].y - Bodies[n].y
                rz = Bodies[x].z - Bodies[n].z

                r3 = (rx**2+ry**2+rz**2)**1.5
                gravconst = G*Bodies[n].M*Bodies[x].M
                fx = -gravconst*rx/r3
                fy = -gravconst*ry/r3
                fz = -gravconst*rz/r3


                # Make a temporary list of the total forces and then add them to get the resulting force
                templistx.append(fx)
                templisty.append(fy)
                templistz.append(fz)

        forcex = sum(templistx)
        forcey = sum(templisty)
        forcez = sum(templistz)
        templistx.clear()
        templisty.clear()
        templistz.clear()

        x = int(Bodies[n].x) + int(Bodies[n].vx) * dt
        y = int(Bodies[n].y) + int(Bodies[n].vx) * dt
        z = int(Bodies[n].z) + int(Bodies[n].vz) * dt

        Bodies[n].x = x
        Bodies[n].y = y
        Bodies[n].z = z

        vx = int(Bodies[n].vx) + forcex/int(Bodies[n].M)*dt
        vy = int(Bodies[n].vy) + forcey/int(Bodies[n].M)*dt
        vz = int(Bodies[n].vz) + forcez/int(Bodies[n].M)*dt

        Bodies[n].vx = vx
        Bodies[n].vy = vy
        Bodies[n].vz = vz

        t += dt




print(Bodies[0].name)
print(Bodies[0].x)
print(Bodies[0].y)
print(Bodies[0].z)


print(Bodies[1].name)
print(Bodies[1].x)
print(Bodies[1].y)
print(Bodies[1].z)

它应该输出类似于此处坐标的内容，但也输出 z 坐标:坐标 1 (41.147123353981485, -2812171.2728945166)坐标 2 (150013715707.77917, 2374319765.821534)

但它输出如下:

Sun 0.0, 0.0, 0.0

Earth 149600000000.0, 0.0, 0.0

注意:问题可能出在 for 循环或数组总和的舍入中，但我不确定。

Answer 1

图片 - 1000 字

您的代码中的直接错误是

• 你计算的力方向错误，应该是 rx = b[n].x-b[x].x 等，或者你需要删除一些行的减号稍后。

• 您在单坐标中的计算会导致复制粘贴错误

  x = int(Bodies[n].x) + int(Bodies[n].vx) * dt
  y = int(Bodies[n].y) + int(Bodies[n].vx) * dt
  z = int(Bodies[n].z) + int(Bodies[n].vz) * dt
  

  在 y 坐标中，您仍然使用 vx。整数值的中间舍入没有意义，它只会稍微降低准确性。

---

我更改了您的代码以使用 numpy 数组作为向量，将加速度计算与欧拉更新分开，在数值模拟期间删除了无意义的整数值舍入，删除了未使用的变量和字段，删除了力的中间变量/加速度计算以直接更新加速度场，更改循环以使用时间来通知一年(或 10 年)已经过去(您的代码以 0.1 天的增量迭代 100 多年，是故意的吗？)，...并添加维纳斯到 body 并添加代码以产生图像，结果见上。

这种螺旋是欧拉方法的典型特征。您可以通过将欧拉更新更改为辛欧拉更新来轻松改进该模式，这意味着首先更新速度并使用新速度计算位置。在其他一切都相同的情况下，这给出了图像

day = 60*60*24
# Constants
G = 6.67408e-11
au = 1.496e11

class CelBody(object):
    # Constants of nature
    # Universal constant of gravitation
    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 (kg)
        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 Bodies

t = 0
dt = 0.1*day

Bodies = [
    CelBody(0, 'Sun', [0, 0, 0], [0, 0, 0], 1.989e30, 'yellow', 10),
    CelBody(1, 'Earth', [-1*au, 0, 0], [0, 29783, 0], 5.9742e24, 'blue', 3),
    CelBody(2, 'Venus', [0, 0.723 * au, 0], [ 35020, 0, 0], 4.8685e24, 'red', 2),
    ]

paths = [ [ b.x[:2].copy() ] for b in Bodies]

# loop over ten astronomical years
v = 0
while t < 10*365.242*day:
    # 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 = sum(rx**2)**1.5
            body.a += -G*other.M*rx/r3

    for n, planet in enumerate(Bodies):
        # use the symplectic Euler method for better conservation of the constants of motion
        planet.v += planet.a*dt
        planet.x += planet.v*dt
        paths[n].append( planet.x[:2].copy() )
        #print("%10s x:%53s v:%53s"%(planet.name,planet.x, planet.v))
    if t > v:
        print("t=%f"%t)
        for b in Bodies: print("%10s %s"%(b.name,b.x))
        v += 30.5*day
    t += dt

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