樹結構在Ｎ體問題中的應用
2008-01-01 00:00:00王武馮仰德遲學斌
計算機應用研究 2008年1期
摘要：N體問題的數值模擬在每個時間步都需要計算每對粒子之間的相互作用，其復雜度為O(N2)。采用樹結構代碼不僅減少了存儲開銷，而且更有利于快速計算和并行劃分。Barnes Hut算法(BHA)和快速多極子方法(FMM)都是基于樹結構的快速算法。BHA可快速計算各點受到的場力，計算復雜度為O(N log N)，但計算精度通常只有1%；FMM通過層次劃分和位勢函數的多極子展開計算各點位勢，其復雜度為O(N)，卻能達到任意精度。數值結果表明，樹結構的并行效果也很好。
關鍵詞：N體問題；樹結構；Barnes Hut算法；快速多極子方法；并行劃分
1N體問題
N體問題實際上是計算場中的多個粒子(如靜電場中的分子、引力場中的星體)之間的相互作用及其運動軌道。以靜電場為例，其基本公式和算法偽代碼描述如下[1，2]:
在恒星動力學中，按照上面的方法模擬N個星體的運動，需要計算各星體所受的力，其復雜度為O(N2)。在重新確定位置后，必須重復以上計算。分子動力學中各粒子的庫侖力計算也是如此。由于恒星系或大分子中粒子的數量很多，需要模擬的時間很長（恒星演化）或時間步很短（粒子運動），數值模擬時特別需要快速算法和并行技術。
2樹結構
許多基于樹結構的算法都是快速算法。下面以二維空間為例，說明四叉樹的結構。三維情況下的八叉樹則完全類似。
先用一個正方形包含所有需要計算的點，再將計算區域四等分。依此類推，直到最小的子區域只包含一個計算點。對計算區域作層次劃分后，其信息存儲在一棵四叉樹中(每個父節點最多有四個子節點)；當點的分布不均勻時，還可以對四叉樹作進一步壓縮簡化。圖1說明了隨機分布于正方形中的十個點作多層組劃分后對應的四叉樹和壓縮四叉樹。
對于更復雜的自適應四叉樹，常用于粒子分布極度不均勻的情形(圖2)。基于樹結構代碼的程序在計算和存儲各個節點信息時會更加方便，而且計算復雜度更低。
1986年Barnes和Hut提出了一種基于樹結構代碼的快速算法(Barnes Hut算法[2])。該算法的總體復雜度為O(N log N)。其計算過程如下：
a)根據粒子分布情況，按照下面的方法構造四叉樹(三維時對應于八叉樹)。理想情況下，即N個點的分布相當均勻，而且最終落在樹的同一層，則構造樹的總體復雜度為O(N log N)；或者根據需要對樹的層數加以限定，即log N≤b，則復雜度為O(b N)。
如何按照點的分布構造四叉樹的簡單算法[2]如下：
procedure QuadtreeBuild
Quadtree={empty}
for i = 1 to n//對所有粒子作循環
QuadInsert(i， root)//將粒子i插入四叉樹相應的位置
end for//四叉樹中可能有很多空葉節點， 但其兄弟節點非空
traverse the tree (via， say， breadth first search)//寬度周游
eliminating empty leaves//去掉空的葉節點
b)通過樹的前序周游(先訪問所有子節點，后訪問父節點，或者稱向上周游)計算每個子區域內部的所有粒子的質心和總質量，存儲在該子區域對應的樹節點上。
c)對每個粒子，通過樹的周游計算它所受到的力。假如在第k層，子區域邊長為D，區域質心和粒子的距離為L。通常限制D/L<θ，常數θ<1控制誤差的大小和計算量[3]。因此BHA只計算某個范圍內的質心對粒子的作用。雖然降低了計算復雜度，卻增加了計算誤差(通常大于1%)。
按照上面的思路， BHA的代碼描述[4]如下：
for(t=0;t { construct Octree (or Quadtree);//構造樹 compute Mass and Center;//計算質量和質心 traverse tree to compute forces;//通過樹的周游計算力 update;//更新各個粒子的信息 } 3快速多極子方法 耶魯大學的L. Greengard和V. Rokhlin在1987年發明了快速多極子方法[1]。它克服了N體問題的瓶頸：將計算復雜度由O(N2)降為O(N)，而且能達到任意精度。因而被美國計算物理學會評為20世紀十大算法之一。FMM的思想是對位勢函數在遠場作多極子展開，再轉換為近場的局部展開。圖3顯示了FMM和BHA的異同。 4結束語 按照常規方法計算N體問題，每個時間步的計算復雜度為O(N2)，限制了求解規模。基于空間層次劃分的樹結構算法，如Barnes Hut算法和FMM，卻將復雜度降低了一個量級。FMM的計算復雜度僅為O(N)，但能得到任意精度，而且通過樹結構的并行劃分得到很高的并行效率。雖然FMM在計算前需作大量的理論分析，但其計算復雜度和數值精度是其他N體問題的計算方法無法比擬的。 參考文獻： ［1］GREENGARD L，ROKHLIN V.A fast algorithm for particle simulations［J］.Journal of Computational Physics，1987，73(2):325－348. [2]BARNES J，HUT P.A hierarchical O(N log N) force calculation algorithm[J].Nature，1986，324(4):446－449. [3]APPEL A W.An efficient program for many body simulations[J].SIAM J Sci Statist Comput，1985，6(1):85 103. [4]WILKINSON B，ALLEN M.并行程序設計[M].陸鑫達，等譯.北京: 機械工業出版社， 2002:103 110. [5]ALURU S，SEVILGEN F.Dynamic compressed hyper octrees with applications to Nbody problem[C]//Proc of Foundations of Software Technology and Theoretical Comp Sci.1999:21－33. [6]HARIHARAN B，ALURU S.Efficient parallel algorithms and software for compressed octrees with applications to hierarchical methods[J].Parallel Computing，2005，31(3－4):311－331. [7]HU Y，JONSSON S L，TENG S H.A data parallel adaptive Nbody method[C]//Proc of the 8th SIAM Conference on Parallel Processing for Scientific Computing.Minneapolis:SIAM Press，1997:19－33. “本文中所涉及到的圖表、注解、公式等內容請以PDF格式閱讀原文”
