天文學家中有這樣一群人，他們既不觀星，也不刷公式，而是通過模擬來研究浩瀚星空。在星團、星系這樣的恒星系統中，往往包含著上百萬顆恒星，規模如此驚人的恒星系統該如何處理？本期“賽先生天文”為你解讀天文學中恒星系統數值模擬的開端與發展。

球狀星團模擬的RGB圖像(Credit: MPA)

撰文 | 李碩(中國科學院國家天文臺)

編輯 | 韓越揚

提起天文學家，我們腦海里最先想到的可能就是一群夜貓子。他們晝伏夜出，沒事喜歡守在望遠鏡前，不知道在搗鼓些啥，總之看上去很厲害就是了。或者，像有些科幻電影里的角色一樣，手里抱本厚厚的《星系動力學》，辦公室黑板上寫滿了連符號都看不懂的公式。當然后者看上去更拉風些。畢竟守著望遠鏡“發呆”的難度系數并不算高，可沒事在黑板上刷公式就不是誰都玩得起的了。

其實，還有那么一小撮天文學家，他們覺得天天盯著望遠鏡和黑板太沒創意，而有了計算機就可以當上帝，想要什么都可以用計算機來模擬。這群喜歡沒事趴在顯示器前，十指噼里啪啦敲個不停的家伙管自己叫：搞數值模擬的。

今天，我們就來看看這個行當是怎么發家的。篇幅有限，我們只聊最“簡單”也是最早開始發展的數值模擬——恒星系統直接動力學模擬中的直接N體數值模擬。

預測恒星的運動

圖1. 著名的球狀星團杜鵑座47，含有數百萬顆恒星(圖源：哈勃拍攝，NASA and Ron Gilliland (Space Telescope Science Institute)；地面拍攝，David Malin, ? Anglo-Australian Observatory)

我們的銀河系是由數千億顆恒星組成的。其中還有大量由眾多恒星聚集而成的星團。既然這些系統都是由恒星構成的，而支配恒星運動的又是萬有引力，那么我們是不是能夠通過計算每顆恒星受到的引力來預測他們的運動呢？

想法挺好，但問題的關鍵是計算量。以最簡單的三體系統為例，要想近似地在紙面上計算，需要把系統演化的過程分成若干間隔盡可能小的時間段(這就是計算機模擬中所謂的時間步長)，然后在每個時間點計算每個天體受到其它天體的引力之和，再根據已有的位置速度“預測”下一時間每個天體的位置速度，循環反復。

聰明的你可能已經意識到了，這類計算的復雜度隨著天體數目是以平方指數增加的。而一個恒星系統，不說像銀河系那樣有幾千億顆恒星，也不說像球狀星團那樣動輒幾十萬顆恒星(不信？可以試著在圖1里面數數)，就算只有十幾顆，計算量也足以讓人吐血了。所以說，沒有計算機，這事也就是想想，真要去做那簡直是瘋了。

“瘋狂”天文學家

圖2. 最早的星系交匯模擬與合并星系NGC2207照片的比較(圖源：Erik Holmberg, 1941, The Astrophysical Journal, 94, 385; NGC2207: ESO)

現在有請“瘋狂”天文學家Erik Holmberg登場。這位仁兄在沒有計算機的年代完成了兩個恒星系統交匯的模擬。圖2左邊是他模擬的恒星系統交匯過程中的形態演化，右圖則是我們拍到的合并星系NGC 2207。

這項工作可是在1941年發表的，還要再等五年，人類第一臺計算機才會誕生。前面提過，沒有計算機幫助這樣的計算可是要出人命的。而Holmberg活到了92歲，顯然并沒在33歲的時候就累死在稿紙堆里。那他是怎么做到的呢？

Holmberg找了74個燈泡來代表兩個星系。對，你沒看錯，74個燈泡！他給不同的燈泡通上了不同的電壓來代表星系中恒星的密度分布。越靠近中心電壓越高，燈泡也就越亮，他通過移動燈泡來模擬星系在交匯過程中的演化。可怎么計算每個燈泡下一步應該挪到哪里去呢？Holmberg找來了光電管測量各處的光強。因為光與引力一樣是遵循距離平方反比衰減的，他巧妙地利用光強的測量代替了引力的計算。

就這樣，史上第一個天體系統動力學演化的模擬就算是完成了。當然，這個模型的分辨率很低——用37個燈泡代替整個星系是有點勉強的。那多加燈泡不行嗎？醒醒吧，剛才我們提到的球狀星團有多少恒星來著？幾十萬顆。幾十萬個燈泡，要開燒烤店么？這畫面太美，實在不敢想。

圖3. Hoerner用來完成首例星團數值模擬的晶體管計算機西門子 2002(圖源：Wikipedia)

計算機誕生后，德國天文學家Sebastian von Hoerner在五十年代末開始了在晶體管計算機上的模擬。當時的計算機沒有顯示器和鍵盤(圖3)。所有程序都要借助穿孔的紙帶輸入計算機，得到的輸出也是一條打滿了孔的紙帶。那時候程序員們每天的工作就是扎眼，代碼不按行來計算而是按捆算。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也就慢慢地退出了歷史舞臺。

圖4. 天龍數值模擬使用的超級計算機Hydra(圖源：馬克思普朗克計算與數據中心)

懸賞一瓶威士忌

九十年后期，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加速時代。

圖5. 天龍數值模擬星團“快照”與球狀星團M13照片(圖源：模擬結果，Wang, L., Monthly Notices of the Royal Astronomical Society, 2016, 458, 1450；M13照片，NASA, ESA, and STScI/AURA)

終于，在2015-2016年，經過與Spurzem、Aarseth等人多年的合作，來自北京大學的博士生王龍完成了Nbody程序的并行GPU版本——Nbody6++GPU，并在先進的通用GPU異構機群上(圖4，可以和圖3比較一下)一舉實現了百萬恒星量級的球狀星團模擬——天龍星團數值模擬，毫無懸念地贏得了Heggie的威士忌。

圖5中可以看到經過可視化處理的數值模擬結果。左圖是模擬結果按星團中不同種類天體呈現的“觀測”效果，中間的圓圖是所有天體疊加的結果，與右圖中Hubble望遠鏡拍攝的球狀星團M13非常接近。而與觀測不同的是，數值模擬可以輕松地區分各種天體，甚至連觀測中基本無法發現的黑洞都能輕松標識。左圖中最右邊的白色背景小圖就是恒星質量級黑洞在星團中的分布。

故事講到這里該告一段落了。但還要多啰嗦兩句，我們只介紹了直接N體方法數值模擬的成長，實際上隨著這一領域的發展，出現了很多方法用以研究恒星系統的演化。即便是N體模擬，在九十年代以后也出現了大量的替代程序，比如Starlab、phi-GPU、Hi-GPU等。篇幅所限就不一一介紹了。

圖6. 銀河系中心恒星軌道，所有恒星都圍繞著在中心不可見的超大質量黑洞(圖源：UCLA Galactic Center Group - W.M. Keck Observatory Laser Team)

最后想說的是，Nbody的故事還遠未到落幕的時刻。我們雖然已經能夠實現星團的模擬，但仍然有很多物理過程無法很好地實現。比如說星系中心的超大質量黑洞(圖6)或球狀星團中心可能存在的中等質量黑洞，現有的N-body程序還無法很好地處理這些極大質量天體與普通天體的相互作用。要想模擬一個有著幾千億顆恒星的星系，我們恐怕還有相當長的路要走，好在我們總能找到幾個喜歡天天坐在電腦前面發呆的家伙來干這事。

致謝：Rainer Spurzem教授為本文的撰寫提供了大量資料，在此深表謝意！

作者簡介

李碩，中國科學院國家天文臺助理研究員。2012年于北京大學物理學院天文學系獲得理學博士學位。研究領域是引力系統演化，集中在超大質量黑洞與星系的共同動力學演化。

與50位技術專家面對面20年技術見證，附贈技術全景圖

總結

以上是生活随笔為你收集整理的matlab模拟三体运动_从灯泡到超级计算机，如何模拟浩瀚星空？| 赛先生的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。