next up previous
Next: 3 FMMの原理 Up: 高速多重極展開法とツリー法---多体シミュレーションのための高速算 法 Fast multipole Previous: 1 はじめに

2 ツリー法の原理

以下、問題を質点系の重力相互作用ということにする。別にクーロン相互作用 やそれ以外の(ポテンシャルでないような)相互作用でも、距離が大き くなれば大きさが小さくなるようなものであれば原理が変わるわけではない。 また、粒子系以外でも、渦糸法による流体計算、偏微分方程式の境界値問題な ど、積分形で書けるものには広く応用できる。そういった研究もなされている。 (例えば[7])

解くべき問題は、以下のような常微分方程式系の数値解を求めるということで ある。

 

ここで Nは粒子の数、 , はそれぞれ粒子 iの位置、質量 であり、 G は重力定数である。

右辺を評価するもっとも単純な方法は、右辺をそのまま計算することである。 粒子数が数百程度であれば、この方法(単に直接計算という)よりも速くする のは難しい。

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

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

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

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

 
図 1: 4分木の構築

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

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

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

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

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

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

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

と書けることになる。

ツリー法を実際に使うには、各節点についてその中の粒子全体の重心の位置と 質量(あるいは多重極展開の展開係数)をあらかじめ求めておく必要がある。 これらはやはり再帰的に定義することができる。つまり、子節点の中の粒子全 体の重心の位置と質量がわかっていれば、親節点の情報は計算できる。多重極 展開の展開係数についても(計算は繁雑であるが原理的には)おなじことであ る。実際の計算プログラムにおいては、再帰手続きで実現してもよいが並列化 などを考えると葉節点から根節点に向かって順に計算していくほうがよい。



Jun Makino
Tue Sep 1 21:46:06 JST 1998