首页 / 技术问答 / Python问答列表 / 问答详情

Python 为什么我的无体模拟器不能打印出3个物体的轨道时间

这是我模拟行星轨道的代码.当我只设置了地球、太阳和木星的天体列表时，我的代码运行良好，并打印出木星和地球轨道的相当准确的时间.然而，当我将土星添加到天体列表中时，我得到的木星和土星轨道的值都是43200秒.然而，奇怪的是，即使在天体列表中有土星，地球的轨道也能准确打印出来.

除此之外，我得到的曲线图是正确的，就像我运行模拟了28年，土星显示没有一个完整的轨道.如果我运行29.5年的模拟，土星just大约完成了它的轨道，这是准确的.所以引力计算仍然在起作用.也许当另一具身体被添加时，指数不知何故出现了故障？但我不确定为什么.

看起来这条线

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()

推荐答案

出现此行为的原因是在这种情况下没有正确触发时间记录的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)

Python相关问答推荐

将numpy数组与空数组相加

除了Python之外，可以替代bare？

Django注释：将时差转换为小数或小数

telegram 机器人API setMyName不起作用

避免循环的最佳方法

用Python获取HTML Span类中的数据

在Arrow上迭代的快速方法.Julia中包含3000万行和25列的表

使用图片生成PDF Django rest框架

pyautogui.locateOnScreen在Linux上的工作方式有所不同

Python panda拆分列保持连续多行

使文本输入中的文本与标签中的文本相同

组/群集按字符串中的子字符串或子字符串中的字符串轮询数据框

实现自定义QWidgets作为QTimeEdit的弹出窗口

移动条情节旁边的半小提琴情节在海运

改进大型数据集的框架性能

在单个对象中解析多个Python数据帧

如何从需要点击/切换的网页中提取表格？

Polars Group by描述扩展

如何将一组组合框重置回无 Select tkinter？

Python—在嵌套列表中添加相同索引的元素，然后计算平均值

实用课程推荐
相关教程推荐
休闲君评测 |  TOOLFK工具网 |  古诗文网 |  良许Linux教程网 |  Python之禅 |  红色石头 |  入党申请书 | 

© 2020 Hi LearnFk 意见反馈 免责声明 关于我们 桂ICP备11002319号-3 友链&广告位+QQ: 1963612630