注意：我已经编辑了原始问题以试图更清楚。

我正在尝试使用 Barnes-Hut 算法模拟引力 nbody 问题，其中物体以随机速度从均匀球体开始。因此，初始角动量几乎为零。我已经用太阳系之类的东西测试了代码，并将其结果与蛮力算法进行了比较，它似乎工作正常。

尽管使用了非常低的时间步长（我使用的是速度 verlet 积分器，如果这很重要），角动量增加了很多，根本不守恒。此图显示了动量的典型增加（稍后我将显示更详细的数据）

您可以看到能量完全守恒，但角动量随着物体的加速而不断增加。所有角动量都是根据集群的实际质心计算的。

当一个物体近距离接触并被驱逐出星团时，角动量和能量都会增加很多：

我尝试了不同的时间步长（从$\delta t= 10^{-8}t_0$到$\delta t=10^{-3}t_0$)和重力软化参数$\epsilon$（从$F_{ij} \propto r_{ij}(\epsilon + r_{ij}^3)^{-1}$） 从$1$光年到$10^{-5}$光年（物体之间的初始分离是$6$ly),并且它们都有相同的一般结果（每行有 8 个系列，但偏差很小）：

线条的突然上升是由于近距离接触而被驱逐的身体引起的。我主要关心的是角动量的初始增加，因为它不会伴随总能量的变化。因为真实的星系有角动量，我想也许我应该从非零角动量开始模拟。

我的问题是：

• 这类问题的标准初始条件是什么？它们涉及一个小的初始角动量，还是使用“旋转”集群开始模拟？

• 在 nbody 模拟中，哪些角动量和能量的变化可以被认为是合理的？

编辑：我没有应用重力软化，但它不会改变结果。

编辑：这是角动量的相关代码，我可以看到一个缺陷，但谁知道......

double acc = 0;
for(int i=0;i<nstars;i++){
    double r[3] = { x[i]-cmx , y[i]-cmy , z[i]-cmz };
    double v[3] = { vx[i] , vy[i] , vz[i] };
    double li[3] = { r[1]*v[2]-r[2]*v[1], -(r[0]*v[2] - r[2]*v[0]),  r[0]*v[1] - r[1]*v[0] };
    // Li = Rᵢ×Vᵢ
    acc += m[i]*norm(li);
    // acc += m⋅R×V
    }
return acc;

cmx, cmy, cmz 是质量的坐标（刚刚检查过，它没有移动一点，它应该在哪里）。使用固定的参考系没有区别。

编辑：刚刚尝试了重力软化（实施得很糟糕）。它有助于避免物体飞走，但对角动量的不断增加没有影响。