python - N体模拟开始时的能量激增
问题描述
我正在尝试使用 n 体模拟来模拟球状星团的坍塌,从所有物体在固定半径球体的随机位置生成的初始速度均为 0 的点开始,我正在调查花费的时间形成一个平衡系统。但是,当粒子最初靠近时运行我的代码时,势能会大幅增加。有没有办法避免这种情况?我正在使用跳跃近似来计算粒子的后续速度和位置。
初始状态
sun_mass = 1.989e30
N = 50
mass = 25 * sun_mass
M = np.full([N],mass)
R = 3.086e+16 # 1pc
epsilon = 0.1 * R
collision_radius = 7e8
V = np.zeros([N,3])
P = np.zeros([N,3])
t = 50000000 * 365 * 24 * 60 * 60
dt = 18000000 * 24 * 60 * 60
print(t/dt)
np.random.seed(54321)
for i in range(N):
#M[i] = np.random.uniform(sun_mass, 100*sun_mass, 1)
phi = np.random.uniform(0,(2*np.pi))
costheta = np.random.uniform(-1,1)
u = np.random.uniform(0,1)
theta = np.arccos( costheta )
r = R * (u) **(1/3)
x = r * np.sin( theta) * np.cos( phi )
y = r * np.sin( theta) * np.sin( phi )
z = r * np.cos( theta )
P[i] = (x,y,z)
运行程序的代码:
def programe(position, mass, velocity, softening, time, dt, R, collision_radius):
#print(len(mass))
no_of_time_steps = (round(time/dt))
#all_positions = np.full((no_of_time_steps, len(mass), 3), 0.0)
all_positions = []
all_velocities = []
#print(all_positions)
#print(len(all_positions[0]))
kinetic_energy = []
potential_energy = []
total_energy = []
#previous_velocity = calc_previous_half_velocity(velocity, calc_acceleration(position, mass, softening), dt)
for i in range(no_of_time_steps):
position, mass, velocity = detect_collisions(position, mass, velocity, collision_radius)
#all_positions[i] = position
all_positions.append(position)
all_velocities.append(velocity)
'graph'
plots = np.round(np.linspace(0,no_of_time_steps,num=500))
for k in range(len(plots)):
if i == plots[k]:
print("test")
#print(i)
graph(position, mass, R, i, dt, k)
'energies'
kinetic_energy.append(calc_kinetic_energy(velocity, mass))
potential_energy.append(calc_potential_energy(position, mass))
total_energy.append(calc_total_energy(position, velocity, mass))
'leap frog'
velocity = calc_next_v_half(position, mass, velocity, softening, dt)
position = calc_next_position(position, mass, velocity, dt)
velocity = calc_next_v_half(position, mass, velocity, softening, dt)
#columns_to_delete = len(all_velocities)
#print(no_of_time_steps)
#print(len(kinetic_energy), len(potential_energy), len(total_energy))
all_positions = np.array(all_positions)
all_velocities = np.array(all_velocities)
#print(all_positions)
graphing(all_positions, all_velocities, mass, time, dt, kinetic_energy, potential_energy, total_energy, no_of_time_steps, R)
#print(len(mass))
return(all_positions, all_velocities, kinetic_energy, potential_energy, total_energy)
蛙跳功能:
'LeapFrog functions'
'Acceleration Calculation'
def calc_acceleration(position, mass, softening):
G = 6.67 * 10**(-11)
N = position.shape[0] # N = number of rows in particle_positions array
acceleration = np.zeros([N,3])
#print(N)
for i in range(N):
#print(i)
for j in range(N):
if i != j:
#print("j", j)
dx = position[i,0] - position[j,0]
dy = position[i,1] - position[j,1]
dz = position[i,2] - position[j,2]
#print(dx,dy,dz)
inv_r3 = ((dx**2 + dy**2 + dz**2 + softening**2)**(-1.5))
acceleration[i,0] += - G * mass[j] * dx * inv_r3
acceleration[i,1] += - G * mass[j] * dy * inv_r3
acceleration[i,2] += - G * mass[j] * dz * inv_r3
return(acceleration) #, print(acceleration))
def calc_previous_half_velocity(velocity, acceleration, dt):
previous_velocity = np.zero_like(velocity)
previous_velocity = velocity - acceleration * dt/2
return(previous_velocity)
def calc_next_v_half(position, mass, velocity, softening, dt):
half_velocity = np.zeros_like(velocity)
half_velocity = velocity + calc_acceleration(position, mass, softening) * dt/2
return(half_velocity)
def calc_next_velocity(position, mass, velocity, softening, dt):
next_velocity = np.zeros_like(velocity)
next_velocity = velocity + calc_acceleration(position, mass, softening) * dt
return(next_velocity)
def calc_next_position(position, mass, velocity, dt):
next_position = np.zeros_like(position)
next_position = position + velocity * dt
return(next_position)
使用的其他功能:
def calc_CoM(position, mass):
sumMR = np.zeros(3)
sumM = 0.0
position = np.array(position)
for i in range(len(mass)):
sumM += mass[i]
sumMR += mass[i] * position[i]
return(sumMR / sumM)
def calc_total_mass2(mass):
total_mass = 0.0
print((mass[0]))
for i in range(len(mass)):
total_mass += mass[i]
return(total_mass)
def calc_total_mass(mass):
total_mass = sum(mass)
return(total_mass)
def calc_CoM_seperation(position, mass, i):
new_mass = np.array(mass)
new_pos = np.array(position)
new_mass[i] = 0.0
new_pos[i] = 0.0
position = np.array(position)
r = np.linalg.norm(calc_CoM(new_pos, new_mass) - position[i])
return(r)
def calc_seperation(p1, p2):
return np.linalg.norm(p1-p2)
解决方案
推荐阅读
- excel - Excel VBA将文件另存为默认文件夹中的word文档
- python - 如何通过数值而不是python中的类频率来排序matplotlib x轴
- python - Python 2.7: Strange Unicode behavior
- python - 使用 for 循环从列表创建过滤列表
- .htaccess - 如何同时使用 gzip 和 http2?
- c# - Visual Studio 2017 社区版即时窗口错误:在运行模式下无法计算表达式
- c++ - C++ 从带有指针数组的函数内部调用函数
- php - 如何清理长 MySQL 查询?
- javascript - 在 dataTables 库中加载数据时添加用于编辑/删除的新列
- javascript - 为什么 contenteditable beforeinput 事件不可取消?