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

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

当粒子彼此非常接近时，能量会增加。如果你采取无限小的步骤，一切都会好起来的。但是如果你采取有限大小的步骤，一个粒子会从 space 中的一个点跳到另一个点。如果它刚好跳到另一个粒子旁边，它的排斥力F_x=delta(x)/r^3就会非常大，这代表它增加了它不应该得到的势能。如果步骤被分解成更小的步骤，粒子就会减慢速度并且不会靠得那么近。

我不知道解决方案。但也许如果你检测到某个步骤的能量增加，那么该步骤可以递归地细分为更小的步骤。

Answer 2

你的物理方程应该是

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

其中 dt 是模拟时间步长。您不会在更新粒子的第二个循环中乘以 dt。这是错误的。

另外你的解法不稳定。您必须选择时间步长以保持稳定性。由于您的系统最终会稳定下来，您可以对其进行试验。

所以建议的操作是：

引入一个新变量dt

将第二个循环更改为

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 个粒子来说看起来不错。

Answer 3

我将力乘以 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

或图表

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

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

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

How come this system of particles keeps increasing in energy?

c

physics