Как эта система частиц продолжает расти в энергии?

Question

У меня есть система частиц. У них есть ускорение, скорость и положение.

Когда частица попадает к стене, ее скорость переворачивается. Когда две частицы приближаются друг к другу, они отталкиваются друг от друга с этой силой:

F=1/r^2 

или

F_x=delta(x)/r^3 
F_y=delta(y)/r^3 

Когда система работала, я чувствовал общую скорость всех частиц увеличиваются. Что странно. Частица должна дать свою энергию другой. Итак, полная энергия системы должна оставаться постоянной.

Кинетическая энергия системы равна

E_k=Sigma v^2 

Я держал контроль полной энергии всей системы и распечатать его с помощью cout и я замечаю, что он продолжает расти. Это противоречит консервативности энергии. Где я делаю ошибку внутри кода?

#include <vector> 
#include <cstdlib> 
#include <cmath> 
#include <iostream> 

constexpr size_t N=1000; 

struct Point 
{ 
    double x, y; 
    double v_x, v_y; 
    double a_x, a_y; 
}; 
Point points[N]; 

void next_frame() 
{ 
    double energy=0.0; 

    // calculate forces 
    for(size_t i = 0; i < N; ++i) 
    { 
     double fx=0.0,fy=0.0; 

     for(size_t j = 0; j < N; ++j) 
     { 
      if(i!=j) 
      { 
       double dx=points[i].x-points[j].x; 
       double dy=points[i].y-points[j].y; 
       double r2=dx*dx+dy*dy; 
       if(r2>0.01 && r2<100.0) // avoid nan and also unnecessary computation 
       { 
        // F=1/r^2 
        double r=sqrt(r2); 
        fx+=dx/(r*r*r); 
        fy+=dy/(r*r*r); 
       } 
      } 
     } 

     points[i].a_x=0.01*fx; 
     points[i].a_y=0.01*fy; 
     energy+=points[i].v_x*points[i].v_x+points[i].v_y*points[i].v_y; 
    } 
    std::cout<<energy<<std::endl; 

    for(size_t i = 0; i < N; ++i) 
    { 
     // integrations 
     points[i].v_x += points[i].a_x; 
     points[i].v_y += points[i].a_y; 
     points[i].x += points[i].v_x; 
     points[i].y += points[i].v_y; 

     // wall 
     if(points[i].x < -50.0) 
      points[i].v_x = +std::abs(points[i].v_x); 
     else if(points[i].x > +50.0) 
      points[i].v_x = -std::abs(points[i].v_x); 

     if(points[i].y < -50.0) 
      points[i].v_y = +std::abs(points[i].v_y); 
     else if(points[i].y > +50.0) 
      points[i].v_y = -std::abs(points[i].v_y); 

    } 

} 

int main(int argc, char **argv) 
{ 
    // initialize particles 
    for(size_t i = 0; i < N; ++i) 
    { 
     Point p; 
     p.x = -50 + ((rand() % 1000)/1000.0)*100.0; 
     p.y = -50 + ((rand() % 1000)/1000.0)*100.0; 
     p.a_x=0.0; 
     p.a_y=0.0; 
     p.v_x=0.001*((rand() % 1000)/1000.0-0.5); 
     p.v_y=0.001*((rand() % 1000)/1000.0-0.5); 
     points[i]=p; 
    } 

    while(1) 
    { 
     next_frame(); 
    } 

    return 0; 
} 

Это энергетический профиль относительно итераций:

---

Пожалуйста, не меняя теги этого вопроса.

Если я задаю этот вопрос на физическом форуме, они скажут мне, что это проблема программирования, а не физика.

Answer 1

Я повторно интерпретировал умножение силы на 0,01 с использованием временного шага dt = 0,01. Тогда используемые скорости на самом деле в 0,01 раза больше реальных скоростей. Чтобы извлечь эту неявную обработку шага по времени, инициализировать скорости с коэффициентом в 100 раз больше,

 p.v_x=0.1*((rand() % 1000)/1000.0-0.5); 
     p.v_y=0.1*((rand() % 1000)/1000.0-0.5); 

и удалить фактор между силой и ускорением

 points[i].a_x=fx; 
     points[i].a_y=fy; 

Затем нанесите временный шаг в процессе интеграции , ([Скорость] Верло симплектична Эйлер с несколькими различных начальными значениями. Из-за случайную инициализацию, это не имеет значения в данном случае.)

 points[i].v_x += points[i].a_x*dt; 
     points[i].v_y += points[i].a_y*dt; 
     points[i].x += points[i].v_x*dt; 
     points[i].y += points[i].v_y*dt; 

---

Чтобы избежать особенностей в гладком и почти физическом изменении пути потенциал использования модифицированного радиуса

r2 = dx*dx + dy*dy + 1e-2; r=sqrt(r2); 

Тогда вы можете удалить условную оценку. Добавьте суммирование потенциальной энергии в этом же цикле

  double r2=dx*dx + dy*dy + 1e-2; 
      // V=1/r, F=1/r^2 
      double r=sqrt(r2); 
      fx+=dx/(r*r*r); 
      fy+=dy/(r*r*r); 
      potential += 1/r; 

и на выходе также объединить кинетическую и потенциальную энергию. С этими изменениями я получаю выход как

kin= 1.70606, pot= 29897.4, tot= 29899.1 
kin= 3.28869, pot= 29895.9, tot= 29899.2 
kin= 7.98328, pot= 29891.3, tot= 29899.2 
kin= 15.4178, pot= 29884.1, tot= 29899.5 
kin= 24.9195, pot= 29875,  tot= 29900 
kin= 35.686,  pot= 29864.9, tot= 29900.6 
kin= 47.0385, pot= 29854.2, tot= 29901.3 
kin= 58.5285, pot= 29843.4, tot= 29901.9 
kin= 69.9214, pot= 29832.6, tot= 29902.5 
kin= 81.1222, pot= 29822,  tot= 29903.1 
kin= 92.1124, pot= 29811.5, tot= 29903.6 
kin= 102.946, pot= 29801.1, tot= 29904 
kin= 113.739, pot= 29790.6, tot= 29904.4 
kin= 124.69,  pot= 29779.9, tot= 29904.6 
kin= 136.055, pot= 29768.8, tot= 29904.9 
kin= 147.937, pot= 29757.3, tot= 29905.2 
kin= 160.059, pot= 29745.7, tot= 29905.7 

или как граф

Как можно видеть, что, в то время как кинетическая энергия неуклонно растет, общая движется энергия только очень медленно. Последний может иметь два источника,

• почти точно сохраняются количество методов симплектической интеграции является модифицированной функцией энергии и
• отражения на границах могут ввести небольшие скачки в этой видоизмененной энергии, что сумма до устойчивого если небольшое изменение полной энергии.

Answer 2

Энергия может увеличиваться, когда частицы очень близки друг к другу. Все будет работать нормально, если вы предпримете бесконечно маленькие шаги. Но если вы делаете шаг с конечным размером, частица переходит из одной точки в пространстве в другую. Если произойдет прыжок прямо рядом с другой частицей, его сила отталкивания F_x=delta(x)/r^3 будет очень большой, что представляет собой увеличение потенциальной энергии, которую оно не должно было получить. Если бы этот шаг был разбит на более мелкие ступени, частицы бы замедлились и не подошли так близко.

Я не знаю о решении. Но, возможно, если вы обнаружите увеличение энергии на шаге, тогда этот шаг можно рекурсивно подразделить на более мелкие шаги.

Answer 3

Ваши физические уравнения должны быть

F = dr/r^3 
a = 0.01 * F 
v += a * dt 
x += v * dt 

Где dt находится на временной шаг моделирования. Вы не умножаетесь на dt во втором цикле, который обновляет частицы. Это не правильно.

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

Так предложил действие:

1. ввести новую переменную dt

2. Измените вторую петлю, чтобы выглядеть

   points[i].v_x += dt * points[i].a_x; 
   points[i].v_y += dt * points[i].a_y; 
   points[i].x += dt * points[i].v_x; 
   points[i].y += dt * points[i].v_y; 
   

Я экспериментировал с небольшим числом частиц в 1D и обнаружил, что dt = 0.001 выглядит отлично для моей группы из 100 частицы.