python 中两体问题的集成

我正在尝试使用 python 和 pygame 创建一个二体 Sim 作为更大项目目标的第一阶段，以在屏幕上显示对象。

我目前的主要问题是，轨道卫星在目标行星周围倾斜时它应该处于稳定的320公里圆形轨道上。

我为四种不同的集成制作了四种不同的功能。 Euler、Leapfrog、Verlet 和 RK4。我预计 Euler 和 Leapfrog 会有一些不准确，但 RK4 或 Verlet 不会。我的微积分知识有限，所以我需要一些额外的眼睛来检查。

它们如下：

def leapfrog_integration(satellite, planet, dt): #most accurate under 5 orbits
    # Update velocity by half-step
    satellite.velocity += 0.5 * satellite.acceleration * dt
    # Update position
    satellite.position += (satellite.velocity / 1000)
    # Calculate new acceleration
    satellite.acceleration = satellite.acceleration_due_to_gravity(planet)
    # Update velocity by another half-step
    satellite.velocity += 0.5 * satellite.acceleration * dt

def euler_integration(satellite, planet, dt):
   # satellite.accleration = satellite.calculate_gravity(planet)
    satellite.acceleration = satellite.acceleration_due_to_gravity(planet)
    satellite.velocity += satellite.acceleration * dt
    satellite.position += (satellite.velocity / 1000) #convert into kilometers
    
def verlet_integration(satellite, planet, dt):
    acc_c = (satellite.acceleration_due_to_gravity(planet) / 1000)#convert to km/s
    satellite.velocity = (satellite.position - satellite.previous_position)
    new_pos = 2 * satellite.position - satellite.previous_position + (acc_c * dt) 
    satellite.previous_position = satellite.position #km
    satellite.position = new_pos #km
    satellite.velocity = (satellite.position - satellite.previous_position)

def rk4_integration(satellite, planet, dt):# need to resolve the conversion to km for position. If i remove the DT from the kx_r then it's the excat same as Verlet and Euler
    def get_acceleration(position, velocity):
        temp_mass = PointMass(position,satellite.mass,satellite.radius,satellite.colour,(velocity), np.zeros_like(satellite.acceleration))
        return temp_mass.acceleration_due_to_gravity(planet)
    
    k1_v = dt * get_acceleration(satellite.position, (satellite.velocity))
    k1_r = (satellite.velocity / 1000)
    
    k2_v = dt * get_acceleration(satellite.position + 0.5 * k1_r, (satellite.velocity) + 0.5 * k1_v)
    k2_r = (satellite.velocity + 0.5 * k1_v) / 1000
    
    k3_v = dt * get_acceleration(satellite.position + 0.5 * k2_r, (satellite.velocity) + 0.5 * k2_v)
    k3_r = (satellite.velocity + 0.5 * k2_v) / 1000
    
    k4_v = dt * get_acceleration(satellite.position + 0.5 * k3_r, (satellite.velocity) + 0.5 * k3_v)
    k4_r = (satellite.velocity + 0.5 * k3_v) / 1000
    
    satellite.position +=(k1_r + 2*k2_r + 2*k3_r + k4_r) / 6 
    satellite.velocity +=(k1_v + 2*k2_v + 2*k3_v + k4_v) / 6

给你点质量的类：

class PointMass:
    def __init__(self, position, mass, radius, colour, velocity, acceleration):
        self.position = np.array(position) #in KM
        self.mass = mass #in Kilograms
        self.radius = radius #in meters
        self.colour = colour
        self.velocity = np.array(velocity) #in m/s
        self.acceleration = np.array(acceleration) #in m/s per second
        self.gForce = None #This is in Newtons
        self.previous_position = self.position - (self.velocity / 1000)  # Initialize previous position for Verlet integration

def acceleration_due_to_gravity(self,other):
        dVector = self.position - other.position # distance vector from self to the other point mass in pixels(km)
        distance_km = np.linalg.norm(dVector)  #Compute the Euclidean distance in km

        distance_m = distance_km * 1000 + other.radius #the distance including the radius to the centre in meters
        unit_vector = (dVector / distance_km) #the unit vector for the direction of the force
        acceleration_magnitude = -Constants.mu_earth / distance_m**2
        return acceleration_magnitude * (unit_vector * 1000) #Return the acceleration vector by multiplying the magnitude with the unit vector(converted to meters)
        #the returned acceleration vector is in m/s

初始条件是：

planet = PointMass(
    position=[600.0,400.0,0.0],
    mass=5.9722e24,
    radius=6.371e6,
    velocity=[0.0,0.0,0.0], #in m/s
    acceleration=[0.0,0.0,0.0], #in m/s^2
    colour=[125,100,100]
)

satellite = PointMass(
    position=[280.0,400.0,0.0],
    mass=100,
    radius=1,
    velocity=[0.0,7718.0,0.0], #need an intial velocity or else it'll jsut fall to the central mass
    acceleration=[1.0,0.0,0.0],#added an non-zero acceleration jsut to make sure there's no issues with the integrations.
    colour=[255,255,255]
)

代码有点混乱，因为我一直在尝试新事物，但还没有真正格式化它。

现在我不是 100% 问题出在集成上，但我想我会从那里开始并继续下去。这很可能是我的重力函数、增量时间应用错误、单位转换导致问题或它在屏幕上的放置方式才是问题。

Github 链接

上面是存储库链接，只是为了将内容保留在其中这里很整洁，因为它可以是任何东西。测试等的屏幕截图

现在可能只是期望Python的行为会失去浮点数的准确性，但我对此打了X，因为这更可能是我糟糕的微积分知识。

经过广泛的测试，我DT 为 0.021，这是增量时间的最佳结果。 Leapfrog 似乎对于低轨道计数或高轨道计数都是最准确的。

Euler 和 Verlet 似乎甚至在 100 个轨道以下，其中 Euler 有点小，但是 RK4 似乎不稳定，因为它慢慢减慢，因此轨道变大和更大。

我在每个集成的不同位置都进行了从米到公里的转换，因为它们无法正常工作，具体取决于我放置的位置。

我必须从某些集成部分中删除 DT，因为如果我把它们留在物体中，它们只会稍微螺旋并落到中心质量上。

这里是对代码的分析，以及对使的 N 体模拟更加准确的步骤：

1. 单位一致性： 的模拟中最大的问题是单位不一致。混合使用米和千米来表示位置和速度，这会导致计算中出现错误。 请始终坚持使用一组单位。 我建议在整个模拟中使用米和秒。

2. 牛顿万有引力定律： 的 acceleration_due_to_gravity 函数看起来基本正确。但是，为了清楚起见，以下是使用米作为单位的版本：

def acceleration_due_to_gravity(self, other):
    dVector = self.position - other.position  # distance vector in meters
    distance = np.linalg.norm(dVector)  # distance in meters
    acceleration_magnitude = -Constants.mu_earth / distance**2
    return acceleration_magnitude * (dVector / distance)  # acceleration vector in m/s^2

1. 集成方法：

2. 蛙跳法： 的蛙跳法实现几乎正确。但是，应该在更新位置 之后 更新速度。

   python def leapfrog_integration(satellite, planet, dt): # Update position satellite.position += satellite.velocity * dt # Calculate new acceleration satellite.acceleration = satellite.acceleration_due_to_gravity(planet) # Update velocity satellite.velocity += satellite.acceleration * dt

3. 欧拉法： 的欧拉法实现是正确的。但是，请记住，欧拉法是一种一阶方法，因此在长时间段内会导致较大的误差。

4. Verlet 积分： 的 Verlet 积分实现不正确。Verlet 积分的正确方程式是：

   python def verlet_integration(satellite, planet, dt): satellite.position += satellite.velocity * dt + 0.5 * satellite.acceleration * dt**2 temp_acc = satellite.acceleration satellite.acceleration = satellite.acceleration_due_to_gravity(planet) satellite.velocity += 0.5 * (temp_acc + satellite.acceleration) * dt

5. RK4 积分： 的 RK4 积分实现看起来很复杂，因为它试图处理单位转换。以下是简化的 RK4 实现，在整个过程中使用米和秒：

   ```python def rk4_integration(satellite, planet, dt): def acceleration(pos): return satellite.acceleration_due_to_gravity(planet)

    k1_r = satellite.velocity * dt
    k1_v = acceleration(satellite.position) * dt
   
    k2_r = (satellite.velocity + 0.5 * k1_v) * dt
    k2_v = acceleration(satellite.position + 0.5 * k1_r) * dt
   
    k3_r = (satellite.velocity + 0.5 * k2_v) * dt
    k3_v = acceleration(satellite.position + 0.5 * k2_r) * dt
   
    k4_r = (satellite.velocity + k3_v) * dt
    k4_v = acceleration(satellite.position + k3_r) * dt
   
    satellite.position += (k1_r + 2*k2_r + 2*k3_r + k4_r) / 6
    satellite.velocity += (k1_v + 2*k2_v + 2*k3_v + k4_v) / 6
   

   ```

6. 初始条件： 的初始条件应该反映希望卫星具有的轨道。对于 320 公里的圆形轨道，可以使用以下公式计算所需的轨道速度：

v = sqrt(G * M / r)

其中： * G 是万有引力常数 (6.67430e-11 m^3 kg^-1 s^-2) * M 是地球的质量 (5.9722e24 kg) * r 是轨道半径（地球半径 + 轨道高度，在本例中为 6.371e6 m + 320000 m）

确保将此速度矢量垂直于卫星到地球中心的位移矢量，以实现圆形轨道。

1. 时间步长： 的时间步长 ( dt ) 会影响模拟的准确性和稳定性。较小的时间步长会产生更准确的结果，但会增加计算时间。对于的模拟，我建议从 1 秒或更短的时间步长开始，并根据需要进行调整。

2. 可视化： 在 pygame 中可视化模拟结果时，需要将卫星的位置从米转换为像素。这可以通过使用缩放因子和将原点从左上角移到屏幕中心来实现。

其他提示：

• 在整个代码中添加注释，以解释的代码在做什么。这将使更容易理解和调试的代码。
• 在每个时间步打印卫星的位置和速度，以帮助调试代码并了解模拟的行为。
• 考虑使用 Python 的调试器来逐步执行的代码并找出任何错误。

通过遵循这些步骤，可以提高 N 体模拟的准确性和稳定性。