next up previous
Next: 4. FMMの原理 Up: 研究展望:高速多重極展開法---粒子法への応用を中心として Previous: 2. 重力多体問題

3. ツリー法の原理

ツリー法や高速多重極展開法の基本的な考えは、遠くの粒子を適当にまとめる というものである。例えば遠くの方の粒子のかたまりをそれらの重心で置き換 えてもいいし、高い精度が必要なら多重極展開も考えられる。遠くにいくに従っ てかたまりを大きくしていけば、直観的には一つの粒子への力を で計算でき、全粒子への力であればその N倍の で計算できる ことになる。

もちろん、このようにまとめられるのは相互作用の性質による。つまり、相互 作用が粒子間距離の関数であり、遠くなればゼロに近付く、また、空間微分も 同様にゼロに近付くという性質があれば、ある計算精度を達成しようと思え ば、遠くはまとめて扱い、近くは細かくみるというやりかたができるわけであ る。

しかし、ここで問題なのはどうやってうまく粒子をかたまりに分けるかとい うことであろう。ある粒子への力を計算するのに適した分け方を見つけるには、 すくなくとも一度は全粒子をみないといけない。従って、その計算量は より小さくはない。粒子毎に別の分け方をしたら、それだけで計算量 が になってしまう。

この問題は、すべての粒子に使える汎用の分け方というのをあらかじめ構成す ることができれば解決する。それをするのがツリー構造による階層的な空間分 割ということになる。

図で説明しよう。3次元の図を書くのは面倒なので、2次元で書くことにする (図1)。まず、平面内の適当な領域に粒子が分布しているとしよう。粒子が 有限個なので、そのすべてを含む正方形を考えることが出来る。すべて含んで いればよくて別に最小である必要はないが、普通はそこそこ小さくとる。 これを、まず4つの小正方形(セル)にわける。で、そのそれぞれをさらに4つに わけるというのを再帰的に繰り返していく。

  
図 1: 4分木の構築。上から順に細かく分割している。

通常の実装では、セルの中の粒子が0個または1個になったところで止める。こ の空間(平面)分割に対応したツリー構造を考えると、根に全体の正方形に対 応する節点があり、その下に4つの小正方形、さらにその下にそのそれぞれの 中の4つ と続いていって、ツリーの葉には粒子1つ1つがくることにな る。この方法では、粒子の数密度が高いところでは細かく、そうでないところ では粗いという意味で、適応的な構造が実現されていることに注意してほしい。 3次元にすれば、正方形4個の代わりに立方体8個となるが、それ以外は全く同 じである。 図2 に、図1 に対応するツリーの構造 を示す。

  
図 2: 4分木構造

どうやってこのツリー構造を使って一つの粒子への重力が計算できるか ということを考えてみる。このためには、ツリーのある節点が表す立方体内の 粒子全体からのある粒子への力を計算する方法を与えればよい。(なお、以下 では、節点に対応する立方体、あるいはそのなかの粒子全体のことも単に節点 ということがある)系全体からの力は、単に根節点からの力として計算できる。 一つの節点からの力は以下のように計算できる:

これは再帰的に定義されている。簡単にプログラムを書くには、再帰呼びだし が使える言語を使えばいい。再帰呼びだしが使えない Fortran 77 のような言 語で書く場合には、アルゴリズムを繰り返しの形に変形すればよい。これにつ いては後で少し触れる。

「十分離れている」かどうかの判定には通常は以下のよ うな基準を使う。

ここで d は粒子と節点の重心の距離、 l は節点に対応する立方体の一辺 の長さである。パラメータは見込み角といわれる量であり、計算精度 と計算量を制御する。高い精度を得たければ を小さくすればいいが、 漸近的には に比例して計算量が増える。計算精度をあげるもう 一つの方法は、重心だけを考えるのでなく多重極展開を作っておくことである。

原点中心で半径 a の球内にある粒子が球の外に作るポテンシャルの多重極 展開係数は

で与えられる。ここで、はそれぞれは粒子 i の質量と極座標で表した位置である。また、 l次の球面調和関数で、ルジャンドル陪関数 を使って

と書ける。この展開係数を使うと、位置r>a)でのポ テンシャルは

と書けることになる。

なお、多重極展開を使う場合には、上のような見込み角を使って判定する方法 はあまり良くない。誤差項を評価して分割するかどうかを判定するのが望まし い。

ツリー法を実際に使うには、各節点についてその中の粒子全体の重心の位置と 質量(あるいは多重極展開の展開係数)をあらかじめ求めておく必要がある。 これらはやはり再帰的に定義することができる。つまり、子節点の中の粒子全 体の重心の位置と質量がわかっていれば、親節点の情報は計算できる。多重極 展開の展開係数についても(計算は繁雑であるが原理的には)おなじことであ る。

ここではとりあえず原理について述べたが、実際にプログラムにする際には計 算量の最適化、ベクトル化、並列化等考えるべきことが多々ある。これらにつ いては後で述べる。



Jun Makino
Thu Jul 8 11:44:39 JST 1999