Блог

Я пытаюсь закодировать код моделирования N-body на python, и мне успешно удалось создать систему с участием Солнца, Земли и Юпитера, как показано ниже, используя метод аппроксимации чехарды.

Однако, когда я пытаюсь расширить один и тот же код для 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)
 

Проблема в том, что симплектические методы обладают своими особыми свойствами только до тех пор, пока системы остаются вдали от каких-либо особенностей. Для гравитационной системы это имеет место, если она иерархична, как в солнечной системе с солнцем, планетами и лунами, где все орбиты имеют низкие эксцентриситеты.

Однако, если вы рассматриваете «звездное скопление» с объектами примерно одинаковой массы, вы не получаете эллипсов Кеплера, и вероятность очень близких столкновений становится довольно высокой. Тем более, что ваше начальное условие нулевой скорости приводит к первоначальному свободному падению всех звезд к общему центру тяжести, что также можно увидеть на вашем детальном изображении.

Из-за потенциальной энергии, падающей в сингулярность, кинетическая энергия увеличивается с уменьшением расстояния, поэтому близкие столкновения равны высокой скорости. При постоянном размере шага, как в методе leapfrog-Verlet, частота дискретизации становится слишком малой, чтобы представить кривую, полностью захватить переход. Экономия энергии грубо нарушается, и высокая скорость сохраняется за пределами близкого контакта, что приводит к нефизическому взрыву системы.