我有一个粒子系统。它们有加速度，速度和位置。
当一个粒子撞到墙上时，它的速度就会翻转。当两个粒子相互靠近时，它们会以这种力相互排斥：
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;
}
这是迭代过程中的能量分布：
请避免更改此问题的标签。
如果我在物理论坛上问这个问题，他们会告诉我这是一个编程问题，而不是物理问题。
最佳答案
我用dt=0.01的时间步长重新解释了力乘以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;
然后在集成过程中应用时间步。（[速度]verlet是初始值略有不同的辛euler。由于是随机初始化，在这种情况下这无关紧要。）
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
或作为图表
我们可以看到，当动能稳定增长时，总能量的移动非常缓慢。后者可能有两个来源，
辛积分方法几乎完全保留的量是一个修正的能量函数，并且
边界处的反射可能会在修正能量中引入小跳跃，如果总能量的变化很小，则总能量会稳定下来。
关于c - 这个粒子系统为何不断增加能量？，我们在Stack Overflow上找到一个类似的问题：https://stackoverflow.com/questions/41300459/