3

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

这是我使用的代码：

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)

它应该输出类似于此处的坐标，但也输出 az 坐标： coordinate 1 (41.147123353981485, -2812171.2728945166) coordinate 2 (150013715707.77917, 2374319765.821534)

但它输出以下内容：

太阳 0.0, 0.0, 0.0

地球 149600000000.0, 0.0, 0.0

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

4

2 回答 2

3

图片 - 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 年，这是故意的吗？），......并将 Venus 添加到正文并添加代码以生成图像，结果见上文。

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

在此处输入图像描述

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()
于 2018-12-20T14:09:49.893 回答
1

我认为您问题的核心是您没有将其视为状态引擎。

想象一下“Bodies”是一个完全不可改变的值，它决定了系统在某个时间点的状态：

bodies_at_time_0 = ((sun, position, velocity, mass), (earth, position, velocity, mass))

你得到下一个状态，如下所示：

bodies_at_time_1 = apply_euler_method_for_one_tick( bodies_at_time_0 )

因此，您的“身体”一次完全固定，您下次计算一个全新的“身体”。在计算中，您总是使用输入中的数据，这就是它们现在所在的位置。您正在做的是移动一些东西，然后根据错误的数字计算将其他东西移动到哪里（因为您已经移动了其他东西）。

一旦你确定你的函数使用输入状态，并返回一个输出状态，你就可以更容易地分解它：

# advance all bodies one time interval, using their frozen state 
def compute(bodies):
    new_bodies = []
    for body in bodies:
        new_bodies.append(compute_one_body(body, bodies))
    return new_bodies

# figure out where one body will move to, return its new state
def compute_one_body(start, bodies):
    end = math stuff using the fixed state in bodies
    return end

# MAIN
bodies = initial_state
for timepoint in whatever:
    bodies = compute(bodies)

我喜欢将元组用于此类事情，以避免意外更改其他范围内的列表（因为列表是可变的）。

于 2018-12-17T15:58:16.943 回答