c - 为什么这个粒子系统的能量不断增加？

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;

然后在积分过程中应用时间步长。 ([Velocity] Verlet 是辛欧拉，初始值略有不同。由于随机初始化，这在这种情况下无关紧要。)

        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

或作为图表

正如您所看到的，虽然动能稳定增长，但总能量的移动非常缓慢。后者可能有两个来源，

• 几乎精确保存的辛积分方法的量是修正的能量函数和
• 边界处的反射可能会在修改后的能量中引入小的跳跃，如果总能量发生微小变化，则总计会达到稳定。