提问人：morgy2190 提问时间：11/15/2023 最后编辑：morgy2190 更新时间：11/15/2023 访问量：70
为什么我的 Nbody 模拟器没有打印出超过 3 个天体的轨道时间？
Why is my Nbody simulator not printing out orbital times past 3 bodies?
问：
这是我模拟行星轨道的代码。当我只用地球、太阳和木星设置了天体列表时，我的代码运行良好，并打印出木星和地球轨道的相当准确的时间。但是，当我将土星添加到天体列表中时，我得到木星和土星轨道的值为 43200 秒。但奇怪的是，即使土星在天体列表中，地球的轨道也能准确地打印出来。
除此之外，我得到的图是正确的，就好像我运行了 28 年的模拟，土星显示没有完整的轨道。如果我运行29.5年的模拟，土星几乎完成了它的轨道，这是准确的。所以引力计算仍然有效。也许当添加另一个主体时，索引以某种方式失败？但我不知道为什么。
似乎是这条线
simulation.run(30*u.year,6*u.hr)
是造成这种情况的原因。每当我将时间步长更改为大约 24*u.hr 时，它都会打印出仅针对木星的准确相邻时间，但不够准确。当我将时间步长更改为 10 天时，它会打印出土星轨道周期的值 10600。然而，这还不够准确。
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
#making a class for Celestial Objects
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
#3d array for position of body in 3d space in AU
self.pos=pos_vec
#3d array for velocity of body in 3d space in km/s
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.cgs.value
else:
return self.mass
v_earth=(((c.G*1.98892E30)/1.495978707E11)**0.5)/1000
v_jupiter=(((c.G*1.98892E30)/7.779089276E11)**0.5)/1000
v_saturn=(((c.G*1.98892E30)/1.421179772E12)**0.5)/1000
#set up first instance of a celestial object, Earth
Earth=CelestialObjects(name='Earth',
pos_vec=np.array([0,1,0])*u.AU,
vel_vec=np.array([v_earth.value,0,0])*u.km/u.s,
mass=1.0*c.M_earth)
#set up second instance of a celestial object, the Sun
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)
Jupiter=CelestialObjects(name='Jupiter',
pos_vec=np.array([0,5.2,0])*u.AU,
vel_vec=np.array([v_jupiter.value,0,0])*u.km/u.s,
mass=317.8*c.M_earth)
Saturn=CelestialObjects(name='Saturn',
pos_vec=np.array([0,9.5,0])*u.AU,
vel_vec=np.array([v_saturn.value,0,0])*u.km/u.s,
mass=95.16*c.M_earth)
bodies=[Earth,Sun,Jupiter,Saturn]
#making a class for system
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()
orbit_completed=False
orbit_start_time=clock_time
initial_position=self.quant_vec[:3]
min_distance=float('inf')
min_distance_time=0
count=0
for step in range(nsteps):
sys.stdout.flush()
sys.stdout.write('Integrating: step = {}/{}| Simulation Time = {}'.format(step,nsteps,round(clock_time,3))+'\r')
y_new=self.rk4(0,dt)
self.history.append(y_new)
self.quant_vec=y_new
clock_time+=dt
#checking if planet has completed an orbit
current_position=self.quant_vec[:3] #**THIS IS WHERE ORBIT TIME CALCULATED** To explain, this is earths position, the first three elements of the vector, the next three elements are its velocity, and then the next three are the suns position vectors and so on. Saturns index would be self.quant_vec[18:21].
distance_to_initial=np.linalg.norm(current_position-initial_position)
if distance_to_initial<min_distance and orbit_completed is False:
min_distance=distance_to_initial
min_distance_time=clock_time
if count==1:
orbit_completed=True
count+=1
runtime=time.time()-start_time
print(clock_time)
print('\n')
print('Simulation completed in {} seconds'.format(runtime))
print(f'Minimum distance reached at time: {min_distance_time:.3f} seconds. Minimum distance: {min_distance:.3e}')
self.history=np.array(self.history)
def nbody_solver(t,y,masses):
N_bodies=int(len(y)/6)
solved_vector=np.zeros(y.size)
distance=[]
for i in range(N_bodies):
ioffset=i * 6
for j in range(N_bodies):
joffset=j * 6
solved_vector[ioffset]=y[ioffset+3]
solved_vector[ioffset+1]=y[ioffset+4]
solved_vector[ioffset+2]=y[ioffset+5]
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*masses[j]/r**3)*dx
ay=(-c.G.cgs*masses[j]/r**3)*dy
az=(-c.G.cgs*masses[j]/r**3)*dz
ax=ax.value
ay=ay.value
az=az.value
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(30*u.year,12*u.hr)
earth_position = simulation.history[:, :3] # Extracting position for Earth
sun_position = simulation.history[:, 6:9] # Extracting position for Sun
jupiter_position=simulation.history[:, 12:15] #Jupiter position
saturn_position=simulation.history[:, 18:21] #Saturn position
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
# Plot the trajectories
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.plot(jupiter_position[:, 0], jupiter_position[:, 1], jupiter_position[:, 2], label='Jupiter')
ax.plot(saturn_position[:, 0], saturn_position[:, 1], saturn_position[:, 2], label='Saturn')
# Add labels and title
ax.set_xlabel('X (AU)')
ax.set_ylabel('Y (AU)')
ax.set_zlabel('Z (AU)')
ax.set_title('Trajectories of Earth and Jupiter Around the Sun')
ax.scatter([0], [0], [0], marker='o', color='yellow', s=200, label='Sun') # Marking the Sun at the origin
ax.scatter(earth_position[0, 0], earth_position[0, 1], earth_position[0, 2], marker='o', color='blue', s=50, label='Earth') # Marking the initial position of Earth
ax.scatter(jupiter_position[0, 0], jupiter_position[0, 1], jupiter_position[0, 2], marker='o', color='green', s=100, label='Jupiter') # Marking the initial position of Earth
ax.scatter(saturn_position[0, 0], saturn_position[0, 1], saturn_position[0, 2], marker='o', color='red', s=80, label='Saturn') # Marking the initial position of Earth
# Add a legend
ax.legend()
# Show the plot
plt.show()
答：
1赞
Christian Karcher
11/15/2023
#1
此行为的原因是在这种情况下未正确触发时间记录的条件。if distance_to_initial<min_distance
造成这种情况的具体原因是“min_distance”，在某些情况下（例如，由于步距增加），它可能会被“跳过”，然后只使用第一个min_distance_time值，而从未“完成轨道”。
为了解决这个问题，我不会使用min_distance，而是检查初始位置是否在前一个位置和当前位置之间。
下面是实现此方法的 run 方法：
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()
initial_position = self.quant_vec[18:21]
min_distance = None
min_distance_time = None
for step in range(nsteps):
sys.stdout.flush()
sys.stdout.write(
"Integrating: step = {}/{}| Simulation Time = {}".format(
step, nsteps, round(clock_time, 3)
)
+ "\r"
)
y_new = self.rk4(0, dt)
y_old = self.quant_vec
self.history.append(y_new)
self.quant_vec = y_new
clock_time += dt
current_position = self.quant_vec[18:21]
last_position = y_old[18:21]
if step == 0 or min_distance_time is not None:
# do not check orbit criteria for initial step
# or if orbit already completed
continue
# checking if planet has completed an orbit, i.e. if origin is between current and last position
last_to_initial = np.linalg.norm(last_position - initial_position)
current_to_initial = np.linalg.norm(current_position - initial_position)
current_to_last = np.linalg.norm(current_position - last_position)
has_passed_initial_position = (
last_to_initial <= current_to_last
and current_to_initial <= current_to_last
)
if has_passed_initial_position:
min_distance = current_to_initial
min_distance_time = clock_time
runtime = time.time() - start_time
print("\n")
print("Simulation completed in {} seconds".format(runtime))
print(
f"Minimum distance reached at time: {min_distance_time:.3f} seconds. Minimum distance: {min_distance:.3e}"
)
self.history = np.array(self.history)
评论
min_distance