Python双体模拟异常：卫星无法维持稳定圆轨道求助

我正在开发基于Python和Pygame的双体物理模拟项目，作为更大项目的第一阶段，目标是在屏幕上展示天体运动。目前核心问题是：本该维持320km稳定圆轨道的卫星，却围绕行星出现振荡现象。

我实现了四种积分算法：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]
)

代码因不断调试略显杂乱。不确定问题是否出在积分算法上，也可能是重力函数、delta time应用错误、单位转换问题或屏幕渲染逻辑导致。

已排除Python浮点数精度问题，更倾向于是自身微积分知识不足导致。经测试：

• 最优delta time为0.021
• Leapfrog算法在低/高轨道次数下精度最高
• Euler和Verlet在100次轨道内表现相近，Euler略有偏差
• RK4则不稳定，卫星逐渐减速导致轨道不断变大
• 不同积分算法中米转千米的转换位置不同，部分算法中若保留dt会导致卫星螺旋下落。

内容的提问来源于stack exchange，提问作者Nzjeux