我正在尝试模拟许多物体之间的引力相互作用。我正在使用直接 PP 力计算和具有可变步长的四阶辛积分器。

模拟的能量由这个函数计算：

real energy(body* system,int n){
  int i,k,j;
  real energy_var=0;
  for(i=0;i<n;i++){
    energy_var += 0.5*system[i].mass*pow(mag(system[i].vel),2);
    for(j=i+1;j<n;j++){
      if(i==j){continue;}
      real position_vect[3]={0,0,0};

      for(k=0;k<3;k++){
    position_vect[k] = system[i].pos[k] - system[j].pos[k];}

      real distance = mag(position_vect);

      energy_var += -Gravitational_constant*system[i].mass*system[j].mass/distance;
    }
  }
  return energy_var;}

在 100 个物体的测试运行中，除了一个异常之外，能量（每 100000 次计算）似乎是恒定的。

Energy = -3.23197e+41
Energy = -3.23197e+41
Energy = -3.23197e+41
Energy = -4.24807e+41
Energy = -3.23197e+41
Energy = -3.23197e+41
Energy = -3.23197e+41
Energy = -3.23197e+41
Energy = -3.23197e+41
Energy = -3.23197e+41

如果这是一个数字或编码错误，即时间步长太高，那么我不希望能量返回到相同的常数值。是否存在可能导致此类错误的已知现象？