假设在时间 t，对于每个粒子，你有：
P position
V speed
以及粒子 A(i) 和 A(j) 之间的 N*(N-1)/2 信息数组，其中 i < j；您使用对称性来评估上三角矩阵而不是完整的 N*(N-1) 网格。
MAT[i][j] = { dx, dy, dz, sx, sy, sz }.
这意味着对于粒子 j，粒子 j 的距离由三个分量 dx、dy 和 dz 组成；以及与 dt 相乘的 delta-vee，即 sx、sy、sz。
要移动到瞬间 t+dt，您可以根据所有粒子的速度暂时更新所有粒子的位置
px[i] += dx[i] // px,py,pz make up vector P; dx,dy,dz is vector V premultiplied by dt
py[i] += dy[i] // Which means that we could use "particle 0" as a fixed origin
pz[i] += dz[i] // except you can't collide with the origin, since it's virtual
然后检查整个 N*(N-1)/2 数组并初步计算每对粒子之间的新相对距离。
dx1 = dx + sx
dy1 = dy + sy
dz1 = dz + sz
DN = dx1*dx1+dy1*dy1+dz1*dz1 # This is the new distance
如果 DN < D^2 且粒子的直径为 D，则您在刚刚过去的 dt 中发生了碰撞。
然后，您可以准确计算发生这种情况的位置，即计算碰撞的确切 d't，您可以根据旧距离平方 D2 (dx*dx+dy*dy+dz*dz) 和新 DN 来计算：
d't = [(SQRT(D2)-D)/(SQRT(D2)-SQRT(DN))]*dt
（以在时间 dt 内覆盖距离 SQRT(D2)-SQRT(DN) 的速度减少从 SQRT(D2) 到 D 的距离所需的时间）。这使得从粒子 i 的参考系中看到的粒子 j 没有“过冲”的假设。
这是一个更繁重的计算，但只有在发生碰撞时才需要它。
知道d't，并且d"t = dt-d't，就可以用dx*d't/dt等对Pi和Pj重复位置计算，得到粒子i和j在瞬间的准确位置P碰撞；您更新速度，然后将其整合到剩余的 d"t 并获得时间 dt 结束时的位置。
请注意，如果我们在此处停止，如果发生三粒子碰撞，此方法将中断，并且仅处理两粒子碰撞。
因此，我们不运行计算，而是标记粒子 (i,j) 在 d't 发生碰撞，并且在运行结束时，我们保存发生碰撞的最小 d't 以及发生在谁之间。
即，假设我们检查粒子 25 和 110 并在 0.7 dt 处发现碰撞；然后我们在 0.3 dt 处发现 110 和 139 之间的碰撞。没有早于 0.3 dt 的碰撞。
我们进入碰撞更新阶段，“碰撞”110 和 139 并更新它们的位置和速度。然后对每个 (i, 110) 和 (i, 139) 重复 2*(N-2) 计算。
我们会发现可能仍然与粒子 25 发生碰撞，但现在在 0.5 dt 处，或者说，在 0.9 dt 处可能会在 139 和 80 之间发生另一次碰撞。0.5 dt 是新的最小值，所以我们在 25 到 110 之间重复碰撞计算，然后重复，每次碰撞的算法都会稍微“减速”。
如此实施后，现在唯一的风险是“鬼碰撞”，即粒子在时间 t-dt 距离目标的直径为 D >，而在时间 t在另一侧的直径为 D > 。
你只能通过选择一个 dt 来避免这种情况，这样任何粒子在任何给定的 dt 中都不会超过其自身直径的一半。实际上，您可以使用基于最快粒子速度的自适应 dt。鬼扫视碰撞仍然是可能的；进一步的改进是根据任意两个粒子之间的最近距离来减少 dt。
这样，算法在碰撞附近确实会大大减慢，但在不太可能发生碰撞时会大大加快速度。如果两个粒子之间的最小距离（我们在循环期间几乎没有成本计算）使得最快的粒子（我们也几乎没有成本发现）不能在不到 50 dts 内覆盖它，那就是 4900 ％的速度增加就在那里。
无论如何，在无碰撞的通用情况下，我们现在已经完成了五次求和（实际上更像是三十四次，由于数组索引），每个粒子对的三个乘积和几个分配。如果我们包括 (k,k) 对来考虑粒子更新本身，那么到目前为止我们有一个很好的成本近似值。
这种方法的优点是 O(N^2) - 它随粒子数量缩放 - 而不是 O(M^3) - 随所涉及的空间体积缩放。
我希望现代处理器上的 C 程序能够实时管理数万个数量级的粒子。
PS：这实际上与 Nicolas Repiquet 的方法非常相似，包括在多次碰撞的 4D 附近减速的必要性。