【并行计算】宇宙天体模拟——Barnes-Hut 算法

并行计算实验之对 N 体问题的O(N*logN)数值模拟。

在卢老师的《并行计算》课上接触到了一个很有趣的问题:N体问题,学习了Barnes-Hut 算法。在完成课内实验的基础上,我利用three.js前端库写了一个效果相当炫酷的星系展示模型:星系展示模型。效果如下:


以下记录一下这次实验中用到的技术。

Barnes-Hut算法

在物理学中,n体问题是预测一组天体在重力(万有引力)作用下相互作用的个体运动的问题。对n体问题的模拟可以帮助我们计算和预测太阳,月亮,行星和恒星的运动。假如我们可以对n体问题进行求解,又能够精确地测量出所有天体的质量、相对位置和速度,理论上我们就可以预测出未来这些天体之间的运动状态,进而对宇宙的进程进行预测。在20世纪,了解球状星团系统的动力学成为一个重要的n体问题。广义相对论中的n体问题 要解决起来要困难得多。

17世纪,Johann Bernoulli 已经将两个物体之间的运动问题,也就是两体问题完全解决。但直到如今,对于三个物体之间的三体问题,依然只能进行数值模拟。著名的小说《三体》就是由此得名的。

虽然完整的n体问题难以得到计算结果,但我们可以通过数值模拟的方式对这个问题进行一些模拟。以三体问题为例,它可以表示为一个十八阶微分方程组:


对这个微分方程组进行数值模拟,就可以计算出天体各个时刻的属性了。

稍微解释一下数值模拟的,我们对时间进行分割,然后通过每个时刻天体中的位置和质量,通过相对论公式和万有引力公式来计算下一时刻各个天体之间的受力关系,然后再通过牛顿第二定律计算下一时刻天体的运动状态。通过这些运动状态,我们可以更新下一时刻天体的位置,然后进行下一轮计算。

在实际的计算中,由于18阶微分方程需要的计算复杂度太高,我们将相对论对系统的影响忽略,将整个过程的方程组简化,此时的方程组如下:


简单地用C++实现了这种数值模拟,并将结果输出。代码如下:


可以看到,用这种方式计算星体的复杂度是O(N^2)的,在10000个星球的时候计算速度已经非常缓慢。在实际的问题中,有时会需要对数十万、上百万个对象进行模拟,这样的复杂度显然无法满足要求。因此引出了下面的Barnes-Hut算法

Barnes-Hut算法

Barnes-Hut算法是专门针对N体问题的模拟过程设计的一个算法,它可以将一次模拟过程的复杂度变为O(N*logN)。

该算法的基本思想是:通过八叉树将整个空间划分为立方体。对于每个星球,由于距离较近的星球对它的影响较大,距离较远的星球对它的影响较小,因此只需要对它附近的星球进行单独处理,然后将较远的星球的信息压缩在树中的虚拟节点中。这样就可以将每个星球的计算量压缩到O(log N)级别。

对Barnes-Hut算法的详细介绍可以参见Wikipedia的页面:https://en.wikipedia.org/wiki/Barnes%E2%80%93Hut_simulation。

以下是对Barnes-Hut算法的实现:


界面展示

为了将计算出的数据比较好地展示出来,也是为了快速地可视化开发,我采用了three.js前端框架来展示三维的画面。相关代码如下:


总结

在这次的实验中用到了许多并行计算、递归与分治的算法知识,在前端展示代码时也学到了很多js和3d绘图的技巧。由于学期末时间匆忙,这篇文章写得比较含糊,等有空一定会多加完善。