为何我的N体求解器中地球沿直线运动而非轨道运行？

问题分析与修正 #

1. 核心问题：初始速度方向错误 #

地球初始位置在(0,1,0)AU（y轴正方向），要进入绕太阳的圆周轨道，速度方向需要垂直于位置矢量（即x轴方向），但你设置的速度是(0,30,0)km/s（y轴方向），这会导致地球沿y轴直线运动，受引力减速后反向，而非做圆周运动。

2. 次要问题：代码逻辑冗余与不规范 #

• 微分方程求解器中，位置导数（速度）的赋值被重复执行多次，逻辑冗余；
• RK4调用时固定传入时间参数0，不符合规范（虽对自治系统无影响，但需修正）；
• 绘图时模拟数据用厘米单位，却标注AU，单位不匹配导致显示异常。

修正后的完整代码 #

import numpy as np 
import matplotlib.pyplot as plt 
import astropy.units as u 
import astropy.constants as c 
import sys 
import time
from mpl_toolkits.mplot3d import Axes3D

# 创建天体类
class CelestialObjects():
    def __init__(self,mass,pos_vec,vel_vec,name=None, has_units=True):
        self.name=name
        self.has_units=has_units
        if self.has_units:
            self.mass=mass.cgs
            self.pos=pos_vec.cgs.value
            self.vel=vel_vec.cgs.value
        else:
            self.mass=mass 
            self.pos=pos_vec 
            self.vel=vel_vec
        
    def return_vec(self):
        return np.concatenate((self.pos,self.vel))
    def return_name(self):
        return self.name
    def return_mass(self):
        if self.has_units:
            return self.mass.value
        else:
            return self.mass
  

# 创建天体实例：修正地球速度方向为x轴（垂直于位置矢量）
Earth=CelestialObjects(name='Earth',
                       pos_vec=np.array([0,1,0])*u.AU,
                       vel_vec=np.array([30,0,0])*u.km/u.s,
                       mass=1.0*c.M_earth)
Sun=CelestialObjects(name='Sun',
                     pos_vec=np.array([0,0,0])*u.AU,
                     vel_vec=np.array([0,0,0])*u.km/u.s,
                     mass=1*u.Msun)
bodies=[Earth,Sun]

# 创建模拟系统类
class Simulation():
    def __init__(self,bodies,has_units=True):
        self.has_units=has_units
        self.bodies=bodies
        self.Nbodies=len(self.bodies)
        self.Ndim=6
        self.quant_vec=np.concatenate(np.array([i.return_vec() for i in self.bodies]))
        self.mass_vec=np.array([i.return_mass() for i in self.bodies])
        self.name_vec=[i.return_name() for i in self.bodies]
        
    def set_diff_eqs(self,calc_diff_eqs,**kwargs):
        self.diff_eqs_kwargs=kwargs
        self.calc_diff_eqs=calc_diff_eqs
    
    def rk4(self,t,dt):
        k1= dt* self.calc_diff_eqs(t,self.quant_vec,self.mass_vec,**self.diff_eqs_kwargs)
        k2=dt*self.calc_diff_eqs(t+dt*0.5,self.quant_vec+0.5*k1,self.mass_vec,**self.diff_eqs_kwargs)
        k3=dt*self.calc_diff_eqs(t+dt*0.5,self.quant_vec+0.5*k2,self.mass_vec,**self.diff_eqs_kwargs)
        k4=dt*self.calc_diff_eqs(t+dt,self.quant_vec+k3,self.mass_vec,**self.diff_eqs_kwargs)
        
        y_new=self.quant_vec+((k1+2*k2+2*k3+k4)/6)
        return y_new
    
    def run(self,T,dt,t0=0):
        if not hasattr(self,'calc_diff_eqs'):
            raise AttributeError('You must set a diff eq solver first.')
        if self.has_units:
            try:
                _=t0.unit
            except:
                t0=(t0*T.unit).cgs.value
            T=T.cgs.value
            dt=dt.cgs.value
        
        self.history=[self.quant_vec]
        clock_time=t0
        nsteps=int((T-t0)/dt)
        start_time=time.time()
        for step in range(nsteps):
            sys.stdout.flush()
            sys.stdout.write('Integrating: step = {}/{}| Simulation Time = {}'.format(step,nsteps,round(clock_time,3))+'\r')
            # 修正：传入当前模拟时间而非固定0
            y_new=self.rk4(clock_time,dt)
            self.history.append(y_new)
            self.quant_vec=y_new
            clock_time+=dt
        runtime=time.time()-start_time
        print('\n')
        print('Simulation completed in {} seconds'.format(runtime))
        self.history=np.array(self.history)
    
def nbody_solver(t,y,masses):
    N_bodies=int(len(y)/6)
    solved_vector=np.zeros(y.size)
    
    for i in range(N_bodies):
        ioffset=i * 6
        # 将位置导数赋值移到j循环外，避免重复执行
        solved_vector[ioffset] = y[ioffset+3]
        solved_vector[ioffset+1] = y[ioffset+4]
        solved_vector[ioffset+2] = y[ioffset+5]
        
        for j in range(N_bodies):
            joffset=j * 6
            if i != j:
                dx= y[ioffset]-y[joffset]
                dy=y[ioffset+1]-y[joffset+1]
                dz=y[ioffset+2]-y[joffset+2]
                r=(dx**2+dy**2+dz**2)**0.5
                ax=(-c.G.cgs.value*masses[j]/r**3)*dx
                ay=(-c.G.cgs.value*masses[j]/r**3)*dy
                az=(-c.G.cgs.value*masses[j]/r**3)*dz
                solved_vector[ioffset+3]+=ax
                solved_vector[ioffset+4]+=ay
                solved_vector[ioffset+5]+=az
    return solved_vector

simulation=Simulation(bodies)
simulation.set_diff_eqs(nbody_solver)
simulation.run(365*u.day,1*u.hr)

# 将厘米单位转换为AU，匹配绘图标注
au_to_cm = (1*u.AU).cgs.value
earth_position = simulation.history[:, :3] / au_to_cm  
sun_position = simulation.history[:, 6:9] / au_to_cm    

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

ax.plot(earth_position[:, 0], earth_position[:, 1], earth_position[:, 2], label='Earth')
ax.plot(sun_position[:, 0], sun_position[:, 1], sun_position[:, 2], label='Sun')

ax.set_xlabel('X (AU)')
ax.set_ylabel('Y (AU)')
ax.set_zlabel('Z (AU)')
ax.set_title('Trajectories of Earth and Sun')
ax.legend()

plt.show()

修正说明 #

1. 速度方向修正：将地球速度改为x轴方向，与位置矢量垂直，满足圆周运动的动力学条件；
2. 代码逻辑优化：把位置导数的赋值移到循环外，消除冗余操作；
3. 规范RK4调用：传入当前模拟时间，符合数值积分的规范流程；
4. 单位统一：绘图前将模拟数据从厘米转换为AU，确保显示与标注一致。

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