我正在尝试在 python 中编写 N 体模拟代码，并成功地使用越级近似方法成功地生成了一个涉及太阳、地球和木星的系统，如下所示。
但是，当我尝试为 N 个质量相同且速度为零的物体扩展相同的代码时，我没有得到系统形成的预期结果。取而代之的是，在最初相互吸引后身体散开的地方会产生以下情况。
无论使用多少初始粒子，都会复制相同的模式。
第二张图片只是第一张图片的放大版，显示他们最初是相互吸引的。
让我相信错误一定在于我的初始条件：
N = 3
mass = 1e30
R = 1e10
V = np.zeros([N,3])
M = np.full([N],mass)
P = np.random.uniform(-R, R, (N,3))
epsilon = 0.1 * R
加速度计算：
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)
蛙跳功能：
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_position(position, mass, velocity, dt):
next_position = np.zeros_like(position)
next_position = position + velocity * dt
return(next_position)
实际程序功能：
def programe(position, mass, velocity, softening, time, dt):
no_of_time_steps = (round(time/dt))
all_positions = np.full((no_of_time_steps, len(mass), 3), 0.0)
all_velocities = []
kinetic_energy = []
potential_energy = []
total_energy = []
for i in range(no_of_time_steps):
all_positions[i] = position
all_velocities.append(velocity)
'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)
return(all_positions, all_velocities, kinetic_energy, potential_energy, total_energy)