假设在时间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中发生了碰撞。
然后，您可以精确计算出发生碰撞的位置，即从旧距离平方D2（dx * dx + dy * dy + dz * dz）和新的DN计算出精确的d't：它是
d't = [(SQRT(D2)-D)/(SQRT(D2)-SQRT(DN))]*dt
（在速度覆盖距离SQRT（D2）-SQRT（DN）的时间dt内，减少距离从SQRT（D2）到D所需的时间）。这使得假设粒子j从粒子i的参考系中看到时，没有“超前”。
这是一个更繁重的计算，但只有在发生碰撞时才需要它。
知道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。
这样，当碰撞不太可能发生时，算法的速度会大大加快，但在碰撞附近时，它的速度会明显变慢。如果两个粒子之间的最小距离（我们在循环中几乎不花费成本计算）使得最快的粒子（我们也几乎不花费成本找到）不能在少于五十个dt内覆盖它，那么速度将提高4900％。
无论如何，在没有碰撞的通用情况下，我们现在已经对每个粒子对进行了五次求和（实际上由于数组索引而更像三十四次），三次乘积和多次赋值。如果我们包括（k，k）对以考虑粒子更新本身，我们就可以得出迄今为止的成本的良好近似值。
这种方法的优点是O（N ^ 2）-它随着粒子数量而扩展-而不是O（M ^ 3）-它随着涉及的空间体积而扩展。
我期望在现代处理器上运行的C程序能够实时管理数万个粒子。
附言：这实际上与Nicolas Repiquet的方法非常相似，包括在多个碰撞的4D邻域减速的必要性。
