N体仿真python

问题描述

我正在尝试用 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,velocity,softening,dt):
    half_velocity = np.zeros_like(velocity)
    half_velocity = velocity + calc_acceleration(position,softening) * dt/2
    return(half_velocity)
       

def calc_next_position(position,dt):
    next_position = np.zeros_like(position)
    
    next_position = position + velocity * dt
    
    return(next_position)

实际程序功能：

def programe(position,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,dt)    
        position = calc_next_position(position,dt)    
        velocity = calc_next_v_half(position,dt)
        

    return(all_positions,all_velocities,kinetic_energy,potential_energy,total_energy)

解决方法

暂无找到可以解决该程序问题的有效方法，小编努力寻找整理中！

如果你已经找到好的解决方法，欢迎将解决方案带上本链接一起发送给小编。

小编邮箱:dio#foxmail.com (将#修改为@）