next up previous
Next: 3 それぞれの方法の歴史と現状 Up: 高速多重極展開法とツリー法 --- O(N) は Previous: 1 はじめに

2 原理

ツリー法と FMMはどちらももともとは質点系の重力相互作用の高速計算のため に提案されたものである。別にクーロン相互作用やそれ以外の(ポテン シャルでないような)相互作用でも、距離が大きくなれば大きさが小さくなる ようなものであれば原理が変わるわけではない。また、粒子系以外でも、境界 要素法や渦糸法による流体計算など、積分形で書けるものには広く応用でき、 そういった研究もなされている。これらについては、例えば IEEE Computational Science and Enginnering の特集号等[4]を御覧 いただくことにして、ここでは以下、問題を質点系の重力相互作用ということ にする。

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

 

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

右辺を評価するもっとも単純な方法は、右辺をそのまま力任せに計算すること である。粒子数が数百程度であれば、この方法(単に直接計算という)よりも 速くするのは難しい。しかし、直接計算では で計算量が増えていく ので、粒子数がある程度大きくなると直接計算では現実的ではなくなる。 計算量を減らすには、重力を及ぼす遠くの粒子を適当にまとめることが考えら れる。さらに、重力を受けるほうも、近くにある粒子は同じような場を感じる わけなので、まとめて計算することが考えられる。重力を及ぼす方は多重極展 開ができる。受けるほうも、同様に場を(球面)調和関数展開して計算する。 こちらは通常局所展開と呼ばれる。

この両方をやるのが FMMであり、多重極展開はするけれども局所展開は使わな いで多重極展開をそれぞれの粒子の位置で評価するのがツリー法である。

 
図 1: ツリー法(上)とFMM(下)

どちらの場合でも、粒子がある空間を分割することで多重極展開や球面調和関 数展開を使えるようにする。もっとも簡単に、空間を規則的な格子で分割した 場合を考えてみると、粒子は、自分と、自分に隣接するセルに入っている粒子 からの力は直接計算し、それ以外のセルからの力は多重極展開で計算すること になる。粒子数を N、 格子のセルの数を n とすれば、計算量は の程度となり、最適な n をとれば となる。

  
図 2: FMMの概念。もっとも濃く塗ったセルに入っている粒子への力を評価する。 自分とその回りの8個のセルからの寄与は下のレベルで計算し、さらにその回 りの27個に多重 極展開を使う。その上のレベルは親セルが面倒を見る。

このやり方で無駄なところは、自分から見て非常に遠くにある粒子も近くにあ るものも全部同じ大きさのセルにわけていることである。計算精度のことを考 えるなら、遠くにいくほどセルが大きくなってよいし、それに対応して遠くか らの力を受けるセルのほうも大きくなるべきである。この状況を図 2に示した。正方形(3次元なら立方体)を再帰的に分割していく ことでツリー構造を作り、あるセルは、自分とは隣接していないが自分の親と は隣接しているセルからの寄与だけを計算する。図では2次元なのでそのよう なセルは27個であるが、3次元では189個になることに注意してほしい。親セル が受けている力は調和関数の展開の中心をシフトして自分の分とあわせる。こ の時、上から順番にやっていけば、それぞれのセルについて一度だけ親からの 分をシフトすればよい。親から伝わってきたものは、その上のレベルからの寄 与をすでに含んでいるからである(図 3)。なお、多重極展開の ほうは、逆に子セルの中心での展開を中心をシフトして足し合わせることで親 セルの展開を構成できる。

  
図 3: 局所展開のシフトと加算

このようにすると、最下位レベルでのセルの数 n を粒子数 N と同じオー ダーにとることができて、計算量が になる。

ツリー法の場合では、遠くにいくほど大きなセルを使うのは同じだが、1つの 粒子がツリーの各レベルと相互作用するので、計算量がレベルの数に比例し、 になる。

ただし、ここで注意しておくべきことは、セル間の相互作用の計算量は球面調 和関数展開の次数 p の 4 乗に比例するのに対し、セルと粒子の間の相互作 用は で計算できることである。

もっとも、FMM の場合には、上に述べたようなセル数の最適化を行なえば計算 量が から に減る。従って、理論上の計算量オーダーと しては FMM のほうが小さくなっている。実際上はどうかということを次節で 考える。



Jun Makino
Thu Nov 26 22:10:48 JST 1998