为什么我对地球绕太阳运行的模拟会产生不准确的力值？

Question

我正在尝试对地球绕太阳运行进行简单的Python模拟，但我不知道为什么我的代码计算出错误的力值。我注意到轨道半径的大小正在变化 - 尽管这应该模拟圆形轨道。我知道这是因为力矢量的大小正在变化，但我不确定为什么会这样。

我正在尝试使用以下公式来计算力值...

我使用 Numpy 模块和离散时间步长来计算未来的力值。我有一个 Planet 类，它具有名为 update_position() 的方法，该方法将 Planet 对象作为输入，并计算传入的 Planet 对象在行星上创建的力。

TIME_STEP = 0.01
G = 6.67430e-11

class Planet:
    def __init__(self, x_position, y_position, x_velocity, y_velocity, mass):
        self.r_position = np.array([x_position, y_position], dtype = 'float64')
        self.velocity = np.array([x_velocity, y_velocity], dtype = 'float64')
        self.mass = mass


    def update_position(self, other):
        force = -1 * G * self.mass * other.mass * (self.r_position - other.r_position) / np.power(np.linalg.norm(self.r_position-other.r_position),3)
        self.velocity += force / self.mass * TIME_STEP
        self.r_position += self.velocity * TIME_STEP

下面，我用真实数据创建了地球和太阳行星对象。然后我反复调用地球上的 update_position(sun) 方法来计算它的新位置。

earth = Planet(x_position = 149_597_870_700, y_position = 0, x_velocity = 0, y_velocity = 29_784.8, mass = 5.972e+24)
sun = Planet(x_position = 0, y_position = 0, x_velocity = 0,y_velocity = 0,mass = 1.9891e+30)

...

while running:
    earth.update_position(sun)

Answer 1

我不认为你的表述有什么特别错误......除了与地球单个轨道的时间相比，时间步长过小。如果 dt=0.01，您将需要执行超过 30 亿次时间步才能完成一个轨道。我怀疑问题在于，与 x 或 y 本身相比，您在每个时间步添加到 x 或 y 的差异是如此之小，以至于这些浮点加法中的错误变得至关重要。

以下代码 - 时间步长为 1 小时（=3600 秒），确认轨道是圆形的。

请注意，要强制形成圆形轨道，请确保您的初始条件给出正确的向心加速度：GMm/r^2 = mv^2/r（您的看起来不错）。

import numpy as np
import matplotlib.pyplot as plt

TIME_STEP = 60.0 * 60       # one hour
G = 6.67430e-11

class Planet:
    def __init__(self, x_position, y_position, x_velocity, y_velocity, mass):
        self.r_position = np.array([x_position, y_position], dtype = 'float64')
        self.velocity = np.array([x_velocity, y_velocity], dtype = 'float64')
        self.mass = mass

    def update_position(self, other):
        force = -1 * G * self.mass * other.mass * (self.r_position - other.r_position) / np.power(np.linalg.norm(self.r_position-other.r_position),3)
        self.velocity += force / self.mass * TIME_STEP
        self.r_position += self.velocity * TIME_STEP

earth = Planet(x_position = 149_597_870_700, y_position = 0, x_velocity = 0, y_velocity = 29_784.8, mass = 5.972e+24)
sun = Planet(x_position = 0, y_position = 0, x_velocity = 0,y_velocity = 0,mass = 1.9891e+30)

num_t = int( 365.25 * 24 * 3600 / TIME_STEP + 0.5 )   # one year, give or take

x, y = [], []
for n in range( num_t ):
# for n in range( 100 ):
    earth.update_position(sun)
    x.append( earth.r_position[0] )
    y.append( earth.r_position[1] )

plt.figure(figsize=(6,6))
plt.plot( x, y )
plt.show()

输出：