Python 中二体问题的积分

Question

我正在尝试使用 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_intergration(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 intergrations.
    colour=[255,255,255]
)

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

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

Github 链接

上面是仓库链接，只是为了保持这里的东西整洁，因为它可以是任何东西。测试等的屏幕截图

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

经过大量测试，我得到的 DT 为 0.021，这是增量时间的最佳结果。无论是低轨道计数还是高轨道计数，Leapfrog 似乎都是最准确的。

Euler 和 Verlet 似乎甚至低于 100 个轨道，其中 Euler 有点偏离，而 RK4 似乎不稳定，因为它缓慢减速，因此轨道变得越来越大。

我在每次集成的不同位置都进行了从米到公里的转换，因为根据我放置的位置，它们无法正常工作。

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

Answer 1

我不知道这是否能解决所有问题，但您可能需要考虑的一件事是在

dVector

中分配

PointMass.acceleration_due_to_gravity

的方式。

您将其定义为

self.position - other.position

。这与你的评论相反；通常，如果您想要一个从

self

到

other

的向量，您会采用

other.position - self.position

（参见 此图）。

这是因为如果您想要一个从

self

到

other

的向量，您会希望能够通过将

other

添加到

vector

来到达

self

。也就是说，您想要：

self

+

vector

=

other

您可以使用减法属性将其重写为

vector

=

other

-

self

。

交换顺序意味着您的加速度矢量将指向您想要的相反方向，这可能会扰乱您的点质量的移动方式。

TL;DR：将 dVector 分配切换为

dVector = other.position - self.position

并查看是否有任何变化。