Python Gravity Simulator 仅绘制直线

Posted 2023-02-21

【中文标题】Python Gravity Simulator 仅绘制直线

【英文标题】：Python Gravity Simulator only plots straight lines

【发布时间】：2022-01-21 12:24:22

【问题描述】：

我正在尝试使用 py 制作一个基本的重力模拟，但由于某种原因，它会将所有线条绘制成直线，当我遇到困难时，我查看了几个示例，但它们都使用类似的方程式/绘制数据的方法所以不知道我哪里出错了

class Body:

    # A class to initialise a body of mass to plot

    def __init__(self, id, mass, coordinates, velocities):
        self.id = id
        self.coordinates = np.array(coordinates, dtype=float)
        self.v = np.array(velocities, dtype=float)
        self.mass = mass
        self.a = np.zeros(3, dtype=float)
        MOTION_LOG.append("name": self.id, "x":[coordinates[0]], "y": [coordinates[1]], "z": [coordinates[2]])


    # Procedure for checking gravity effects on body
    def gravity(self):
        self.a = 0
        for body in bodies:
            if body != self:
                dist = body.coordinates - self.coordinates
                r = np.sqrt(np.sum(dist**2))
                self.a += (SETTINGS['G'] * body.mass * dist / r**3) ** 2


# Procedure to plot the new coordinates of the body         
    def move(self):
        self.v += self.a * SETTINGS["deltaT"]
        self.coordinates += self.v * SETTINGS['deltaT']

然后我真的模拟了它

# For loop to run a simulation for a specific time set
for step in range(int(SETTINGS["tLimit"] / SETTINGS["deltaT"])):
    SETTINGS['elapsedT'] += SETTINGS['deltaT']
    if SETTINGS['elapsedT'] % SETTINGS["frequency"] == 0:
        prog = ((SETTINGS['elapsedT'] / SETTINGS['tLimit'])*100)//1
        for index, location in enumerate(MOTION_LOG):
            location["x"].append(bodies[index].coordinates[0])
            location["y"].append(bodies[index].coordinates[1])
            location["z"].append(bodies[index].coordinates[2])
        print(f"prog%")

    for body in bodies:
        body.gravity()
    for body in bodies:
        body.move()

然后终于在我做的图表上绘制它

fig = plt.figure()
ax = fig.add_subplot(1,1,1, projection='3d')

for body in MOTION_LOG:
    ax.plot(body["x"], body["y"], body["z"])

plt.show()

不确定是我刚刚做错了加速还是我只是画错了点，但我见过的其他例子并没有太大的不同

示例数据

SETTINGS = 
    'G' : 6.67e-11,
    'deltaT' : 172800,
    'elapsedT' : 0,
    'tLimit' : 315360000,
    "frequency": 1,


MOTION_LOG = []
bodies = []
set_bodies = 
    "Sun": 
        "id": "Sun",
        "mass": 1e20,
        "coordinates": [0, 0, 0],
        "velocities": [0, 0, 0]
    ,
    "Earth": 
        "id": "Earth",
        "mass": 1e3,
        "coordinates": [1e2, -1e2, 0],
        "velocities": [0, 0, 0]
    

for body, body_data in set_bodies.items():
    bodies.append(Body(**body_data))

【问题讨论】：

它应该如何绘制线条而不是直线？

@mkrieger1 我希望以曲线/轨道而不是直线绘制线条，例如来自here

@AJBiffl 啊忘了把它放在嵌入中，现在添加它

我认为这是我的加速度方程错误，每当我检查数值时，加速度似乎非常小，所以身体几乎不会改变位置

【参考方案1】：

首先，你的加速度是平方的。你的行是：

self.a += (SETTINGS['G'] * body.mass * dist / r**3) ** 2

你想要的是：

self.a += (SETTINGS['G'] * body.mass * dist / r**3) 

其次，你的号码完全不正常。永远不要，永远只是猜测数字的好组合用于这样的模拟。要么使用完全自然的单位（重力常数、太阳质量和 AU 都等于 1），要么使用所有实际数字。当你不这样做时，你只是在猜测什么时间步是合适的，而且很难做到这一点。例如，根据您用于太阳质量、重力常数和地球轨道半径的数字，一个时间单位约为 7.7 年。要获得看起来不错的东西，那么您需要大约7e-4 的时间步长并运行到大约1.3 的最后时间，并对地球使用[5000, 5000, 0] 的初始速度（您还需要更改您如何添加到MOTION_LOG，因为小型浮点数不能很好地与 mod (%) 配合使用。但不要那样做——使用合理的数字。当您不必花费数小时来回转换成您实际想要的单位时，您最终会更快乐。

第三，你的“地球”物体没有初始速度。无论您如何修正方程式或单位，它都会落入太阳。

以下是一些使用 SI 单位的大约正确的“实数”数字：

set_bodies = 
    "Sun": 
        "id": "Sun",
        "mass": 2e30,
        "coordinates": [0, 0, 0],
        "velocities": [0, 0, 0]
    ,
    "Earth": 
        "id": "Earth",
        "mass": 6e24,
        "coordinates": [1.5e11, 0, 0],
        "velocities": [0, 3e4, 0]
    

【讨论】：

还有一个简短的说明 - 我通常建议在几乎所有模拟中使用自然单位，尤其是轨道系统，因为您遇到机器精度或数字溢出/下溢问题的风险要小得多，因为它们是都接近于一（而不是天文数字）

无论如何我都不倾向于使用我提供的测试数据，我使用了来自here 的数据，这些数据是我抓取的，但即使我将值更改为更自然的值，我仍然会得到直线图表

你试过我在回答中给出的数字吗？

哦，我忘了提到我发现的另一个错误。请参阅答案中的编辑。

哦，是的，当我粘贴正方形时忘记删除它，我稍微改变了我的等式并记得有人提到要平方它，我必须在将它添加到嵌入之前忘记删除它。我的实际方程式是self.a += -(SETTINGS["G"] * body.mass) / r * dist 我只是无法编辑它。即使有这些数字，我也会得到两条相互排斥的直线