天文学家中有这样一群人，他们既不观星，也不刷公式，而是通过模拟来研究浩瀚星空。在星团、星系这样的恒星系统中，往往包含着上百万颗恒星，规模如此惊人的恒星系统该如何处理？本期“赛先生天文”为你解读天文学中恒星系统数值模拟的开端与发展。
提起天文学家，我们脑海里最先想到的可能就是一群夜猫子。他们昼伏夜出，没事喜欢守在望远镜前，不知道在捣鼓些啥，总之看上去很厉害就是了。
或者，像有些科幻电影里的角色一样，手里抱本厚厚的《星系动力学》，办公室黑板上写满了连符号都看不懂的公式。当然后者看上去更拉风些。其实，还有那么一小撮天文学家，他们觉得天天盯着望远镜和黑板太没创意，而有了计算机就可以当上帝，想要什么都可以用计算机来模拟。这群喜欢没事趴在显示器前，十指噼里啪啦敲个不停的家伙管自己叫：搞数值模拟的。
今天，我们就来看看这个行当是怎么发家的。
篇幅有限，我们只聊最“简单”也是最早开始发展的数值模拟——恒星系统直接动力学模拟中的直接N体数值模拟。我们的银河系是由数千亿颗恒星组成的。其中还有大量由众多恒星聚集而成的星团。既然这些系统都是由恒星构成的，而支配恒星运动的又是万有引力，那么我们是不是能够通过计算每颗恒星受到的引力来预测他们的运动呢？想法挺好，但问题的关键是计算量。
以最简单的三体系统为例，要想近似地在纸面上计算，需要把系统演化的过程分成若干间隔尽可能小的时间段（这就是计算机模拟中所谓的时间步长），然后在每个时间点计算每个天体受到其它天体的引力之和，再根据已有的位置速度“预测”下一时间每个天体的位置速度，循环反复。
聪明的你可能已经意识到了，这类计算的复杂度随着天体数目是以平方指数增加的。而一个恒星系统，不说像银河系那样有几千亿颗恒星，也不说像球状星团那样动辄几十万颗恒星（不信？可以试着在图1里面数数），就算只有十几颗，计算量也足以让人吐血了。所以说，没有计算机，这事也就是想想，真要去做那简直是疯了。
现在有请“疯狂”天文学家Erik Holmberg登场。
这位仁兄在没有计算机的年代完成了两个恒星系统交汇的模拟。这项工作可是在1941年发表的，还要再等五年，人类第一台计算机才会诞生。前面提过，没有计算机帮助这样的计算可是要出人命的。而Holmberg活到了92岁，显然并没在33岁的时候就累死在稿纸堆里。那他是怎么做到的呢？Holmberg找了74个灯泡来代表两个星系。对，你没看错，74个灯泡！他给不同的灯泡通上了不同的电压来代表星系中恒星的密度分布。
越靠近中心电压越高，灯泡也就越亮，他通过移动灯泡来模拟星系在交汇过程中的演化。可怎么计算每个灯泡下一步应该挪到哪里去呢？Holmberg找来了光电管测量各处的光强。因为光与引力一样是遵循距离平方反比衰减的，他巧妙地利用光强的测量代替了引力的计算。就这样，史上第一个天体系统动力学演化的模拟就算是完成了。
计算机诞生后，德国天文学家Sebastian von Hoerner在五十年代末开始了在晶体管计算机上的模拟。当时的计算机没有显示器和键盘。所有程序都要借助穿孔的纸带输入计算机，得到的输出也是一条打满了孔的纸带。那时候程序员们每天的工作就是扎眼，代码不按行来计算而是按捆算。Hoerner在1959年就是用这样的计算机实现了16颗星的数值模拟。梦想总是要有的，万一实现了呢？
Hoerner同学估算了一下，发现按照当时的计算速度，“只要”花1400亿年就能模拟一个球状星团了。
几乎是同时，另一位重要的玩家参与到了这项游戏中，这便是挪威人Sverre Aarseth。这位奇人开创了Nbody系列模拟程序、获得过国际象棋通讯赛国际大师称号、炒股挣钱买过法拉利、50岁时迷上登山并因此冻掉过好几根脚趾、80岁时玩野外漂流，为了骗取工作人员许可谎称自己才60。
最有个性的是，在剑桥大学工作了几十年却一直是普通研究员，相当于现在的博士后。要说学术贡献不够吧，有一颗小行星可是以他命名的，用来表彰其贡献的。他的人生经历简直可以写一本书了。
言归正传，Aarseth在A. Schlüter的建议下为每个天体赋予了独有的计算步长。这就避免了两个密近天体因为步长太小而拖慢整个计算的情况。这一方法在1963年提出，于1986年才最终完善为等级阻塞时间步长方案。
从80年代初开始，Aarseth用FORTRAN语言编写了著名的Nbody程序，并在之后与合作者一道维护更新至今。为了更真实地模拟恒星系统，他们又把目标对准了双星。银河系中双星非常普遍。而由于其轨道过于靠近，需要非常小的时间步长才能保证计算精度。这就导致包含双星的计算会非常缓慢。
借助数学的帮助，Aarseth等人在六七十年代通过一种叫做KS（Kustaanheimo & Stiefel）正则化的坐标变换，将双星的轨道积分转换成了一种简单的谐振子计算。后来Seppo Mikkola与Aarseth又在1992年将这一工具应用在了多个天体上，即所谓的KS链。
同时，为了进一步提高精度，日本天文学家Junichiro Makino与Aarseth在91-92年引入了厄米积分方法，通过泰勒展开的数学变换巧妙地用一阶导数取代了模拟中所需的二三阶导数计算，从而在不增加计算复杂性的前提下有效地提高了计算精度。随着计算机硬件的飞速发展，1996年，德国天文学家Rainer Spurzem与Aarseth一起实现了10000个天体的数值模拟。
而在同一年，Makino借助专门设计的计算设备GRAPE-4实现了32768个天体的模拟。GRAPE是一种专门用来进行引力计算的特殊硬件，很容易进行并行化计算。在GPU异构加速计算兴起之前，GRAPE是最高效的引力系统模拟硬件。但是，GRAPE专精于引力计算，很难被大规模应用于其它领域，因此制造成本高昂。而GPU在这方面拥有与GRAPE类似的能力，且成本更低。
在GPU计算兴起后，GRAPE也就慢慢地退出了历史舞台。
九十年后期，Aarseth的Nbody模拟程序已经发展为Nbody6并加入了各种真实的物理过程。随着高性能计算机群的发展，人们开始考虑将模拟工作交给多个CPU并行完成。1999年，Spurzem在Nbody6的基础上开发出了第一个并行版本Nbody6++。几年之后，一些简单的动力学模拟（包括GRAPE）已经可以处理百万天体量级的计算。
当然这些计算还无法兼顾恒星演化等星团中重要的物理过程。为此，英国数学与天文学家Douglas C. Heggie专门悬赏，将为第一个完成真实球状星团模拟的工作赠送一瓶他珍藏的威士忌。2000年后，GPU加速计算迅速兴起。Spurzem尝试借助德国基金会帮助建造GPU机群未果。他得到的答复是：“我们已经有非常强大的超级计算机，为什么还要GPU机群呢？
”而中国此时为了鼓励GPU计算正准备资助建造若干个GPU超级计算机。就这样Spurzem来到了中国，开始将Nbody++与GPU加速相结合。
与此同时，Aarseth与日本天文学家Keigo Nitadori合作，为Nbody6加入了GPU计算支持，使模拟工作可以在一台有GPU加速卡的工作站上完成。Nbody数值模拟由此进入了GPU加速时代。
终于，在2015-2016年，经过与Spurzem、Aarseth等人多年的合作，来自北京大学的博士生王龙完成了Nbody程序的并行GPU版本——Nbody6++GPU，并在先进的通用GPU异构机群上（图4，可以和图3比较一下）一举实现了百万恒星量级的球状星团模拟——天龙星团数值模拟，毫无悬念地赢得了Heggie的威士忌。
故事讲到这里该告一段落了。
但还要多啰嗦两句，我们只介绍了直接N体方法数值模拟的成长，实际上随着这一领域的发展，出现了很多方法用以研究恒星系统的演化。即便是N体模拟，在九十年代以后也出现了大量的替代程序，比如Starlab、phi-GPU、Hi-GPU等。篇幅所限就不一一介绍了。最后想说的是，Nbody的故事还远未到落幕的时刻。我们虽然已经能够实现星团的模拟，但仍然有很多物理过程无法很好地实现。
比如说星系中心的超大质量黑洞（图6）或球状星团中心可能存在的中等质量黑洞，现有的N-body程序还无法很好地处理这些极大质量天体与普通天体的相互作用。要想模拟一个有着几千亿颗恒星的星系，我们恐怕还有相当长的路要走，好在我们总能找到几个喜欢天天坐在电脑前面发呆的家伙来干这事。