今までは、微分方程式の右辺は与えられている、つまり基本的に全粒子からの 力を合計するとして話をしてきたが、ここでは右辺の全粒子からの力を近似を 使って速く計算するという話をする。
歴史的にも、また現在使われているものとしても、右辺を速く計算するための 方法はたくさんある。が、基本的な考え方としては以下の3つにわけることが できよう。
第一は、各粒子への力を直接計算するのではなく、空間格子を切って格子上で ポアソン方程式を解き、それから補間などで各粒子への力を求める方法である。 境界条件や粒子分布によってどのような格子をとるかは違ってくるが、3次元 の問題では直方体領域で普通に等間隔格子を切ることが多い。これは、ポアソ ン方程式を解くのに高速フーリエ変換を使えるので、計算が速いからである。 普通 PM (particle mesh)法というときにはこのFFTを使う方法を指す。
この方法の利点は速度にあるが、欠点はフーリエ変換を使うために規則的な等 間隔格子でないと使えないというところにある。このために空間分解能に制限 がつき、格子間隔の数倍程度以下では粒子間重力は実効的にソフトになるし、 また格子の方向による非等方性も出てくる。このために、例えば宇宙初期で構 造が一様に近いといった、自己重力系としては限られた状況でしか使えない。
この分解能の問題を改善するアプローチももちろんある。昔からあるのは
(particle-particle particle-mesh)法というもので、近くの粒
子からの寄与は直接計算しようというものである。この方法では、 
ポ
テンシャルを2つに分ける。一つは、
の極限では 
だ
が遠くでもっと速く落ちる近距離力的な振舞いをするもの、もう一つはこの近
距離力を引いた残りである。で、前者は直接計算し、後者は FFT を使って計
算する。
この方法は、プログラミングが比較的容易でありまた精度もあげられる。もっ とも、精度をあげようとすると、遠距離力が格子の非等方性の影響を受けない程度に近距 離力の作用範囲を大きくする必要があり、比較的急速に計算量が増える。また、 格子が等間隔であるために粒子の密度分布が非一様であると計算量が非常に急 激に増大する。
これに対し、いわゆる AMR (adaptive mesh refinment) を使うという話も最 近はある。これは、密度が高いとかポテンシャルが急に変わるとかいう ところでは格子を細かくしようというものである。AMR にしてしまうと普通に はFFTが使えないので、 CG法等の反復解法でポアソン方程式を解くことになる。
AMR と 
を組み合わせて近くの粒子からの力も正確に計算しよう
という話もあることはある。[CTP95]但し、 AMR の場合にはメッ
シュを使って計算した力の 
ポテンシャルからのずれを正確に評価する
のは困難になるので、理論的にはむつかしい問題がある。
第二の方法は、基本的には第一の方法と同じようにポアソン方程式を解くが、 格子を切って差分法で解くのではなくなんらかの直交関数系に展開する。といっ ても、FFTでもそのような展開をしているわけで区別が明確とはいい難いが、 ここでは実際に考えるのは球面調和関数展開である。
銀河とか星団とかは、おおむね丸い形をしているので、それを表すのに球面調 和関数展開を使うのは考え方としては非常に自然である。とはいえ、実際にこ の方法を使うに当たってはいろいろ問題がある。根本的な問題は計算量である。 FFT のようなうまい性質がないので、普通にやると計算量は展開項数と粒子数 の積に比例する。従って、角度方向や半径方向の分解能を粒子間隔程度にとる と、結局展開項数は粒子数の程度になり計算量は粒子数の2乗程度になる。従っ て、項数をあまり多くとることは現実的ではない。
実際的な問題としては、 r 方向の展開をどうするかということがある。r 方向の積分について直交する関数ならなんでもいいわけではあるが、現実問題 としては上の項数の問題があるので、あまりなるべく少ない項数で系を表現し たい。このために、もともとの初期条件の分布に近い式を 0次にもってきて直 交関数系を作るような怪しい方法がいくつか提案されている。 [HO92]しかし、項数を制限するということは、系の可 能な進化がその有限項で表されるもののみに制限されるということでもあり、 こういうことをやった時に本当に追跡したい進化が表現できているかどうかは 明らかではなくなる。
もう一つの方法は、半径方向は展開しないということである。つまり、ある半 径での重力ポテンシャルを、内側からの寄与と外側からの寄与に分けると、そ れぞれは多重極展開で表現できる。内側からの寄与は中心に近い粒子から順 番に計算できるし、外側からの寄与はもっとも外の粒子から順に計算できる。
計算量の観点からは、この方法は計算量が角度方向の項数と粒子数の積ですみ、 しかも半径方向については展開による計算誤差が入らないので、3方向に展開 する方法よりも明らかに優れている。このために、1980年代には単独の銀河の 大粒子数計算にはよく使われた。[Vil82,AW86]但し、並列計算を考えると、順番に計算しな いといけないというのが嬉しくないので最近はあまり使われていないようであ る。といっても、プロセッサ数がよほど多くなければ並列化できないところは たいしたことはないので、使われなくなった本当の理由は良くわからない。 この方法の問題は、単独で球に近い銀河でも、やはりAMR と同じくどのよ うな誤差が入っているのか良くわからないということである。
と、いうことで、以下の話の本題であるツリー法系の計算法の話にはいる。ツ リー法の基本的な考え方は、計算量を減らすために重力を及ぼす遠くの粒子を 適当にまとめるということである。まとめるのに、階層的な木構造を使うこと で、どの粒子に対しても小さなコストでうまいまとめ方を構成する。
さらに、重力を受けるほうも、近くにある粒子は同じような場を感じる わけなので、まとめて計算することが考えられる。重力を及ぼす方は多重極展 開ができる。受けるほうも、同様に場を（球面）調和関数展開して計算する。 こちらは通常局所展開と呼ばれる。
この両方をやるのが FMMであり、多重極展開はするけれども局所展開は使わな いで多重極展開をそれぞれの粒子の位置で評価するのがツリー法である。
どちらの場合でも、粒子がある空間を分割することで多重極展開や球面調和関 数展開を使えるようにする。同じ球面調和関数展開を使うにしても、系が「丸 い」ということを使って一箇所で展開するのではなく、分割したものをそれぞ れ展開するので系が丸くなくても使えるわけである。
以下、具体的にどのようにツリー構造を構成し、それによってどのように相互 作用を計算するのかを簡単にまとめる。
図で説明しよう。3次元の図は書きにくいし見る方も見にくいであろうから、2 次元で書くことにする（図10）。まず、平面内の適当な領域に粒子が分布して いるとしよう。粒子が有限個なので、そのすべてを含む正方形を考えることが 出来る。すべて含んでいればよいので、別に最小である必要はないが、普通は そこそこ小さくとる。これを、まず4つの小正方形(セル)にわける。で、その それぞれをさらに4つにわけるというのを再帰的に繰り返していく。
粒子が0個または1個しかなくなったところで止める。この
空間（平面）分割に対応したツリー構造を考えると、根に全体の正方形に対応
する節点があり、その下に4つの小正方形、さらにその下にそのそれぞれの中
の4つ
と続いていって、ツリーの葉には粒子1つ1つがくることになる。
この方法では、粒子の数密度が高いところでは細かく、そうでないところでは
粗いという意味で、適応的な構造が実現されていることに注意してほしい。
3次元にすれば、正方形4個の代わりに立方体8個となるが、それ以外は全く同
じである。
まず、どうやってこのツリー構造を使って一つの粒子への重力が計算できるか ということを考えてみる。このためには、ツリーのある節点が表す立方体内の 粒子全体からのある粒子への力を計算する方法を与えればよい。（なお、以下 では、節点に対応する立方体、あるいはそのなかの粒子全体のことも単に節点 ということがある）系全体からの力は、単に根節点からの力として計算できる。 一つの節点からの力を以下のように計算することにする
「十分離れている」かどうかの判定には通常は以下のよ うな基準を使う。
ここで d は粒子と節点の重心の距離、 l は節点に対応する立方体の一辺
の長さである。パラメータ
は見込み角といわれる量であり、計算精度
と計算量を制御する。高い精度を得たければ 
を小さくすればいいが、
漸近的には 
に比例して計算量が増える。計算精度をあげるもう
一つの方法は、重心だけを考えるのでなく多重極展開を作っておくことである。
原点中心で半径 a の球内にある粒子が球の外に作るポテンシャルの多重極 展開係数は
で与えられる。ここで、
と
はそれぞれは粒子 i の質量と極座標で表した位置である。また、
は l次の球面調和関数で、ルジャンドル陪関数
を使って
と書ける。この展開係数を使うと、位置
（r>a）でのポ
テンシャルは
と書けることになる。
ツリー法を実際に使うには、各節点についてその中の粒子全体の重心の位置と 質量（あるいは多重極展開の展開係数）をあらかじめ求めておく必要がある。 これらはやはり再帰的に定義することができる。つまり、子節点の中の粒子全 体の重心の位置と質量がわかっていれば、親節点の情報は計算できる。多重極 展開の展開係数についても（計算は繁雑であるが原理的には）おなじことであ る。実際の計算プログラムにおいては、再帰手続きで実現してもよいが並列化 などを考えると葉節点から根節点に向かって順に計算していくほうがよい。
ツリー法の原理から FMM の原理にはほんの一歩である。ツリー法では、ある
粒子への力を計算するのに、遠くの粒子からの力はまとめて計算することにし
た。相互作用は対称的であるので、逆にある粒子から遠くの粒子への力をまと
めて計算するということも考えられる。両方を同時にやれば、ある一群の粒子
から別の一群の粒子への力をまとめて計算することになる(図9)。このように、
受ける側もまとめてやることによって、以下に説明するようにツリー法の
から 
へ計算量のオーダーを下げることができる。この
ためにこの方法には「高速」多重極展開法という名前がついている。
力を受けるほうは、重力ポテンシャルのテイラー展開を求めておいて、それを 各粒子の位置で評価するという操作が必要になる。もちろん、実際にはポテン シャルは調和関数になっているので、球面調和関数で展開することで項数をテ イラー展開に比べて減らすことができる。この展開のことを、以下局所展開と 呼ぶ。
この方法は、Rokhlin と Greengard によってまず2次元の場合に提案され
[GR87]、すぐに3次元に拡張された
[GR88]。考え方は2次元でも3次元でも同じであるので、
例によって図は2次元で示す（図11）。ツリー法では適応的なツリー構造を扱っ
たが、 FMM ではとりあえず一様なツリー構造を考える。すなわち、ツリーの
レベルを k とすると、元の正方形を 
個の小正方形に分割する。なお、
適応的なツリー構造の場合のアルゴリズムは例えば[BHER95]を
みられたい。ツリー法の場合と同様に、あらかじめ各セルに含まれる粒子（一
つも無いこともありえるが）が作るポテンシャルの多重極展開は準備しておく。
図 11: FMMの概念。もっとも濃く塗ったセルに入っている粒子への力を評価する。
自分とその回りの8個のセルからの寄与は直接計算し、さらにその回りの27個は多重
極展開を使う。その上のレベルは親セルが面倒を見る。
FMMでは、一つの粒子への力を以下のように分割する
ここで、 
は直接計算する分で、自分のいるセルと、それに隣接した
セルの粒子からの力である。これらには多重極展開は使えないので、直接計算
することになる。残るのが 
であるが、これはさらにツリーのレベル毎
に別に計算する。つまり、あるセル i に対し、(i)そのセルと同じレベルに
あって、(ii)その親セルとその隣接8セルの合計9セルの子であって、(iii)問
題のセルには隣接していない、27 (=36- 9)個のセルを考える。この27個の
セルに含まれる粒子が作る重力ポテンシャルの局所展開の和を、セル iの中
心で評価しておく。トップレベル（系全体）と、そのすぐ下のレベルでは、親
のレベルでは隣接していて自分のレベルでは隣接していないものは存在しない
ので、計算の必要はない。
このようにして各レベルで重力場を計算しておくと、あとは粒子の位置で、そ れぞれのレベルでの局所展開を評価して合計すればいい。が、別に粒子毎に全レベ ルでの展開を評価しなくても、上のレベルの展開は一段下に移してやって（展 開の中心をシフトして）まとめておけば、粒子は全部まとめた局所展開を一度 だけ評価すればいい。この、局所展開をシフトしながら下に落していくのは、 最初に多重極展開をやはりシフトしながら上にあげていくのとちょうど対称的 な操作になる。
図12のように、4レベルのセルの中心での展開を合計する場合を考えてみる。ま ず、一番上のレベルのセルの中心での局所展開を、4つの子セルの中心に展開 中心をシフトしてそれぞれのレベルでの重力場の局所展開と足し合わせる。さ らに、それぞれの子セルで、合計した局所展開をまたそれらの4つの子セルの 中心にシフトし、またそこでの局所展開と足し合わせる。これを最下層のセル に達するまで繰り返すと、各セルの中心での、自分と自分に隣接するセルに入っ ている粒子を除いた他のすべての粒子による力の局所展開がもとまる。粒子一 つへの力は、この最下層での局所展開を粒子の位置で評価したものと、隣接セ ルの粒子からの力の合計ということになる。
ここで計算量のオーダーを考えてみると、すべての操作の計算量が粒子数ま
たはセルの数にしか比例していないことがわかる。セル数は粒子数に比例する
ので、結局計算量のオーダーが 
ということになるわけで
ある。
FMM とツリー法のどちらも、まだ10年少々の歴史しかない新しい方法である。 しかし、どちらも急速に研究が進み、応用分野も粒子系のほか、境界要素法、 あるいは変わったところでは非常に高次の球面調和関数展開や、分散メモリの 計算機上での FFT など、さまざまな分野に広まりつつある。基本的な考え方 は、遠くとの相互作用はまとめるという単純なものなので、重力・クーロン相 互作用に限らず非常にさまざまな応用がありえるわけである。基本的には、相 互作用が重ね合わせ可能でかつ遠距離に届き、といっても遠くにいくに従って 弱くなるようなものであればなんでも FMM やツリーの考え方は使える。
ここでは、特に筆者の専門である重力多体系の観点から過去の歴史を簡単に振 り返ってみたい。
ツリー法を初めて実現したのは、プリンストン大学の学部生 A. Appel であっ
た[App85]。これは1981年の学部の卒業論文であったということだが、
1985 年になるまで公表論文にはならなかった。彼の
方法は、ツリーを最近接粒子をまとめていくことによってボトムアップに構成
していくものであり、現在広く使われている8分木をつかうものではなかった。
8分木を使うツリー法を提案したのは当時プリンストン高等研究所にいた
J. E. Barnes と P. Hut である[BH86]。この方法に基づいて、
その後の数年間にベクトル型スーパーコンピュータへの効率的な実装について
の研究と分散メモリ型計算機への実装の研究が精力的に進められ、銀河や銀河
団、あるいは宇宙の大規模構造などの形成、進化のシミュレーションに広く使
われるようになった[WS92]。特に分散メモリの超並列計算機で
は、 
を越えるような粒子数のシミュレーションが可能になってきてい
る。
このように広く使われることになった理由の一つは、実装が比較的簡単であっ たことである。ツリー法の実装のほとんどは、4重極モーメントまでしか扱 わない。このため、多重極展開や球面調和関数といってもそれほど面倒な 計算式が出てくるわけではない。高い計算精度が必要であれば、離れたセルも ある程度まで分割して相互作用を計算する。もちろん、非常に高い計算精度が 必要になれば、より高次のモーメントをとり入れた方が計算量が減る。しかし、 多くの問題で、粒子数が有限であることからくる2体緩和 で計算精度がきまっているので、相互作用の計算精度をあげるよりも、計算精 度を下げても粒子数を増やしたほうがよりよい結果が得られるのである。この ために、かなり広い応用で、4重極モーメントよりも精度をあげることには意 味がない。これについては後でもう少し詳しく議論する。
このように、理論的にも実験的にも複雑なことをして精度をあげる必要がなかっ たために、ベクトル化や並列化、特に適応的な木構造の実現やロードバランシ ングなどに努力が傾注されてきた。これらの努力により高速で実用的なプログ ラムが利用可能になってきている。
FMM はツリー法とほぼ同時期にイェール大学の Greengard と Rokhlin により
提案された[GR87,GR88]。これは多重極展開と局所展
開をまともに使うものであった。その後、計算法の改良について非常に多様な
提案がなされてきた。例えば多重極展開から局所展開への変換については、
FFTを使って高速化する方法[EJ96]や、一旦座標系を回転し、シ
フトの方向に z 軸をもってくることで計算量を 
に減らす方法
[WHG96]などが提案されている。この他にも最近 Greengard
らによって提案された新しい方法[GR97]もある。これらを
使うことで、理論的な計算コストは 
よりも若干小さくなっている。
しかし、このような多様な提案があるということは、逆にいえばまだ決定版と
いうべきものはなく、多数の改善案が提案されているがそれらのうち一つが他
のものよりも圧倒的によいというわけではないということでもある。さらに、
このような方法は一般に高い精度を要求する場合にしか改善につながらない。
実際のところは、天文シミュレーションにおいては現在のところ FMM は全く 使われていない。このように、理論上の計算量はともかく、現実の問題として は FMM があまり使われていない一つの理由は FMM の実装が面倒であるという ことである。ツリーでは木構造を作った後は多重極モーメント（といってもほ とんどの場合せいぜい 4 重極）を計算し、あとはそれを使って各粒子への力 を計算する再帰的な関数を書けばおしまいだが、 FMM の場合には多重極展開 を局所展開に変換するややこしい式と、局所展開をシフトする式をプログラム しないといけない。しかも、特に多重極展開を局所展開に変換するところはま ともにやったのでは計算量が多いし、といって減らそうと思うとまたいろいろ 面倒なことをしないといけない。
もちろん、それだけやって実際に同じ計算精度で速くなるとか同じ計算時間で
より精度がよくなるとかいうことがあれば、自分で書かないまでも誰かが作っ
たプログラムを使おうかという気にもなるが、現在のところ実際に精度と計算
速度を比較した研究では、少なくとも比較的計算精度が低いところではツリー
の方が速いという結果が得られている。例えば、適応的な木構造の場合にツリー
法とFMMの計算速度を比較した研究[BS97]では、比較的精度が
低い（ポテンシャルの相対精度が 
程度）場合、FMMはツリー法に比
べて有意に遅いという結果になっている。となると難しいことをわざわざしな
いほうが得であろう。
もっとも、現在のところ、ツリー法に比べて FMM の結果があまりよくないの は少なくとも部分的には FMM の実装があまりよくないからであると理論的に は思われる。例えば、前節の説明では、多重極展開を局所展開に変換するとこ ろであるセルは自分と同じ大きさのセルからの力を受ける。で、どちらも同じ 次数で展開する。誤差の上限の評価ではこれは同じオーダーで係数も同程度だ が、実際に入ってくる誤差を考えると必ずしもそうはならない。というのは、 多重極展開のほうは、あるセルに入っている沢山の粒子からの力の合計である ので、高次の項は主にセル内の密度の揺らぎから来ている。このために、中の 粒子数が大きければ揺らぎが小さくなり、その結果誤差も小さくなる。ところ が、局所展開ではあくまでもセル内の一箇所での評価なので、展開の打ちきり 誤差がそのまま入ってくる。中に沢山粒子があっても、それぞれに独立に誤差 がでるだけである。このために、多重極展開と局所展開で実際の誤差に対する 寄与が同程度であるためには、局所展開のほうの次数をかなり高くするか、あ るいは局所展開のほうのセルを小さくする必要があると考えられる。ところが、 現在のところ大抵の実装ではこういった効果まで考慮して計算量を最適化しよ うとはしていない。
なお、2次元の場合には3次元に比べて FMM の計算量は大きく減り、実用性の
高い方法となる。これは、多重極展開の項数が 
から p に減るだけで
なく、幾何学の違いのためにセル間相互作用の数が減るためである。3次元の
場合は一つのセルが 
個のセルと相互作用するが、2次元の場合
は
個であり、7倍計算量が違うことになる。なお、もちろんツリー
法でも同じ理屈で減る。
幾何学の違いは、空間が3(2)次元であっても実際に計算する領域 は2(1)次元である境界要素法の場合にも有利に働くことを注意しておきたい。 このために、FMMやツリー法は境界要素法には特に有効である。
FMM の計算コストを減少させる方法についてはいくつか簡単に述べたが、ここ では実装を簡単にする（必ずしも計算量は減らない）方法を2つ紹介しておく。 これらは、プログラムが簡単に実現できるという意味で、 FMM が広く使われ るためには重要なものと考えられる。
Anderson [And92]は球面調和関数を使わないで FMM を実現する 画期的に安直な方法を提案した。彼のアイディアは、球面調和関数の展開係数 をデータとして使う代わりに、球面でのポテンシャルの値そのものを使おうと いうものである。半径 a の球の中の粒子が外に及ぼすポテンシャルは、球 の表面でのポテンシャルの値がわかっていればラプラス方程式の境界値問題を 解けば与えられる。解はポアソンの公式を使って以下のように書ける
ここで 
は n次ルジャンドル関数である。この積分を、球面上で
の適当な標本点を使った数値積分に置き換える。2次元であれば、円を等分割
すればいいが、3次元ではうまく点をとってその上で球面調和関数の直交性が
保存されるようにする必要がある。これは良く研究されていて、文献[McL63,HS96]を見れば
数値的にそのような点が求められているので、それを使えばよい。
子セルの外の多重極展開をその一つ上のセルでのそれに変換するには、上の式
に従って親セルを含む球面上でのポテンシャルを求めればいい。また、局所展
開も、
での展開が rでの展開に変わるだけで上と同じ式である。さら
に、多重極展開を局所展開に変換するには、上の式をつかって局所展開す
る球面上でポテンシャルを評価すればよい。つまり、上の式だけで、 FMM
に必要な数学がすべて表現されていることになる。
ここで、著者らが最近提案した
法 (pseudo-particle
multipole method)[Mak99,KM01] にも簡単に触れさせていただきたい。基本的な考え方はア
ンダーソンの方法とほとんど同じものであるが、違いは球面上でのポテンシャ
ル分布を使うのではなく球面上の質量分布によって多重極展開を表現するとい
うものである。定式化はアンダーソンの方法とほとんど同じであり、最終的に
は
という形で実粒子の分布から仮想粒子の質量が表現される。ここで
,
はそれぞれ実粒子、仮想粒子の質量であり、 
は展開の中心から見て実粒子と仮想粒子がなす角度、
は実粒子の原点か
らの距離、 r は仮想粒子を置く球の半径、Kは仮想粒子の数である。
この方法は、 FMM よりもむしろツリー法との組み合わせで有効であり、粒子 間相互作用を計算するだけで高次精度のツリー法が実現できる。
もっとも、この方法は計算量という面ではあまり有利とはいい難い。 これは、最小限必要な疑似粒子の数に比べて実際に使う疑似粒子数がかなり多 くなるからである。例えば、もっとも低次の重心（1次のモーメント）だけを 使う場合を考える。疑似粒子とかそういうことを考えなくても、重心に一つ粒 子をおけばいいというのは自明である。ところが、球面調和関数展開してから 1次の項をというようなことを考えると、4個必要になることになる。これが4 重極になると、必要な粒子数は 12 に増える。通常の計算法では、4重極モー メントを計算するのに必要な計算量は粒子一つからの力を計算する場合のせいぜ い2倍程度なので、12というのは普通に考えると現実的な計算量とはいい難い。
必要な疑似粒子数が多いのは、粒子の空間配置を固定しているからである。空 間配置を固定して、質量だけを変えることで、線形計算で疑似粒子を決められ るというメリットはあるものの、粒子あたり4つの自由度のうち一つしか使っ ていないので、最大4倍の粒子が必要になってしまう。
1次の場合には上に述べたように「重心に粒子を置く」というのでいいわけだ が、2次以上に場合にもなんとかならないかというのはそういうわけで検討に 値するであろう。今のところ、2次(4重極)の場合にはそういうものが できるということがわかっている。
2次のモーメントまで考えた時には、独立な係数は9個である。従って、粒子2 個では不可能だが3個ならできそうな気がする。3個の粒子の質量の合計はもち ろん実粒子のそれに等しい必要があるし、重心も一致しなければいけない。こ れで自由度は4個使うので、残っているのは8個である。この8個で四重極モー メントテンソルの独立成分5個を表せればいい。
4重極モーメントテンソルが対角化される方向に座標系をとると、座標系で3 自由度使ったので残りは第一軸方向のモーメントと第二軸方向のモーメント の2つだけである。例えば3粒子の質量を等しくし、それを第一軸か第二軸のど ちらかに底辺が平行な二等辺三角形配置にすることで、対角化の条件は満たさ れるのであとは辺の長さをモーメントの値が等しくなるように決めればいいこ とになる。もっと違う配置ももちろんあり得るが、とりあえず四重極モーメン トを合わせるにはこれで十分といえる。
この方法は実装ができていて、速度はもちろん12個も粒子を使うのに比べれば 断然速い。
4重極（2次）が理論上可能な最小の粒子数で表現できるなら8重極（3次）も、、、 というのが人情だが、今のところそういう方法は発明されていない。これは不 可能とかそういうことではなく、単にまだ誰もそんなことが可能かどうか真面 目に考えたことがないというだけである。
ツリー法の計算精度に関しては、良くわかっていない問題がいくつかある。こ こで問題なのは、例によって「計算された重力に対してある要求した精度を実 現するためにはどうすればよいか」ではなく、「どのような精度になっていれ ばシミュレーション結果は正しいといえるのか」ということである。
例えば、式(133)を使った時に、粒子分布と
を与えれば求まった重力の誤差とかその分布は実験的にも理論的に
もわかる。が、ここで問題なのはではどのような精度が必要かということであ
る。もうちょっというと、それ以前にそもそも重力の誤差というのがシミュレーションの
「正確さ」を測るのに適当な物差しか？という問題が、
あるわけである。
ある粒子の軌道で、ある時刻 tにおいて加速度 
を計算してある誤差
があったとしよう。今、時間積分は無限に正確で誤差がないとし、また
数値積分した結果の軌道のずれは十分小さくて、誤差の非線形性は無視できる
とする。この時、例えばある時間経過した後の速度の誤差はもちろん
を積分し
たものであり、位置の誤差は速度の誤差を積分したものである。
位置や速度自体を正確に計算することにどれほど意味があるかというのは前に
長々と議論したが、まあ例えばエネルギーの誤差とかいうことにしても話は同
じで数値計算法や要求する計算精度によって変わるのは 
であり、さま
ざまな量の誤差は 
の時間積分である。
さて、こんなことはあまりにも当たり前なことではあるが、上で見たように誤
差を決めるのは
の積分であり 
そのものではない。従って、計
算精度を問題にするなら
のある瞬間での値とかその2乗平均値ではなく、
の自己相関関数なりパワースペクトルをみないといけないわけである。
仮に
の瞬間瞬間での値は大きかったとしても、それが短い時間スケール
でランダムに変われば時間積分した時の効果は小さい。逆に瞬間瞬間での値は
小さくても、 長いタイムスケールでしか変わらなければ時間積分した時の効
果は大きいかもしれない。従って、ある瞬間での 
の値自体に意味があ
るとは限らない。
ここで問題なのは、それでは一体どのような尺度を使うべきかということであ るが、これは問題に依存することは明らかであろう。
例えば、ビリアル平衡にある恒星系で、速度分散によって構造が維持されてい るものを考える。この時にツリー法で力を計算すると、ある粒子からの誤差の 相関時間は、その粒子が相互作用計算の時に使われるセルを横断する時間の程 度になる。もちろん、球面調和関数の次数が大きければそれに応じて相関時間 は短くなるが、オーダーとしてはそうなっている。従って、2粒子間の相互作 用の誤差の相関時間は距離にほぼ比例するものと考えられる。というのは、2 粒子間の相対速度は距離に無関係だからである。
あ、ここで暗黙の仮定として、パラメータ 
が一定であるというのを
使っているが、いろいろ最適化をはかった場合でも
の距離依存性は非
常に小さいのでこれはまあ妥当な仮定である。
さて、相関時間が距離に比例するとすると、誤差自体への寄与は距離が近いもの ほど大きくなっても構わないことになる。もちろん、ここにはさらに仮定が必要で、相関時間が 実際に有限であるためには誤差の平均値が 0 でないといけないわけだが、こ れは一応成り立っていると仮定することにする。
さらに、かなり大胆な仮定であるがある粒子への他の2つの粒子からの力の誤 差は無相関であるとする。これは粒子数が大きい極限では明らかに正しくない が、現在可能な計算で実現されている粒子数ではこう仮定して実際上問題はな いようである。
この時には、各粒子からの力の誤差の時間積分の合計をとると、1次の項は 0 になり 2 次の項は結局2体緩和によるエネルギー変化に 2 粒子間の力の相対 誤差程度の大きさの係数がついたものが出てくることになる。 [MIE90,HHM93]従って、粒子の軌道の誤差は2体 緩和によるものよりもずっと小さい。さらに重要なことは、これは結局2体緩 和と同じ距離依存性を持つということなので、近くと遠くの計算誤差が同程度 になっているという、計算量と計算精度の面から見ると最適な状況が実現され ているということである。つまり、粒子分布が一様に近く、粒子のランダム速 度が大きいという状況に関する限り、単純な幾何学的条件が最善であるという ことになる。
ここで、単純な幾何学的条件は賢くなさそうだとか考えて、ある瞬間での力の 誤差を最小にするようなやり方を考えたとしよう。上の議論から明らかなよう に、幾何学的条件を使った場合にはある瞬間での力の誤差に寄与するのは主に 近傍の粒子である。このために、そのような「改善された」やり方では必ず幾 何学的条件に比べて近くをより正確に、遠くをよりいい加減に計算することで 同じ計算量で精度をあげようとすることになる。これで確かに力の誤差は減る が、時間積分した軌道では遠くをいい加減に計算しているのが効いてかえって 悪くなってしまうのである。
もちろん、現実の状況はもっと複雑である。例えば、2つの粒子からの力が誤 差が無相関であるという仮定は空間構造が一様でなければ明らかに正しくない。 一般には小さいスケールでは一様とみなせても大きいスケールではそうではな いわけなので、遠いからといって大きなセルにするのはあまりよろしくないか もしれない。
といっても、多くの銀河のように中心にいくほど密度が高く、遠くまで距離の 冪乗で密度が下がるものを扱う時には、中心から非常に遠いところではそもそ もほとんど質量がないので、相関時間とかいうこと以前に力自体が小さい。こ ういった状況では幾何学的条件は明らかに無駄といえる。
別の例を考える。よく行なわれる宇宙論的な構造形成計算では、宇宙初期には 密度は一様、ランダム速度はほとんどゼロであり、誤差の相関時間は構造形成 が進んで宇宙膨張から切り離された構造ができるまでの間ということになる。 この場合には、相関時間が距離によらないので瞬間での力の誤差を最小にする ような計算法がもっともよい。ところが、同じ宇宙論的計算でも、構造が発達 して銀河や銀河団が出来てくると、時間積分した誤差を最小にするためには瞬 間での力の誤差を最小にするような方法ではいけないことになる。
というようなわけで、現在のところ一つの条件でどんな場合でも使えるといっ た万能な計算精度の決め方はない。結局のところ、計算精度を変えてみて見た い物理現象に影響がなければ多分大丈夫であろうというような実用的な基準で チェックする必要がある。
重力多体系のシミュレーションにおいては、より大きな粒子数を使った計算を するということが極めて重要である。これは、多くの場合に計算精度が2体緩 和で決まっていて、計算精度を改善するには粒子数を増やす以外に効果的な方 法がほとんどないことによる。例えば銀河形成のシミュレーションで2体緩和 によって出来た銀河の構造が変わってしまってはなにを計算しているかわから ないわけだが、銀河の中心付近の高密度な領域、あるいは銀河円盤のように恒 星がおもに回転運動していてランダム速度が小さいような領域では、現在数値計算 で扱える粒子数ではまだ緩和時間が短かく、物理的に意味がある結果を得るこ とは困難である。このような領域では、物理的に意味がある結果を得るために 第一にするべきことはなんらかの方法で扱える粒子数を増やすということにな る。
扱える粒子数は、計算機の速度と計算法の効率で決まるので、計算法の効率を あげるというのと少なくとも同程度には計算機の速度をあげるというのも重要 である。もっとも、後で簡単に触れるような自分で計算機を作るというような ある意味極端なアプローチをとらない限り、計算機は大学なり研究所の計算セ ンターの計算機を時間単位で利用するか、あるいは自分の研究室のワークステー ションやパソコンを使うという話になってそれ以上あまり考える余地がないよ うな気もする。が、これは必ずしもそうでもなくて、実際には使える計算機に 合わせて計算法自体を変更する、あるいは新しいものを開発するといった必要 が出てくる。
特に、ここ10年ほどの高速計算機の発展の方向を見てみると、高いピーク性能 をもつ計算機が安価に供給されるようになったことは事実であるがその上で利 用可能なプログラミング環境はほとんど進歩しておらず、というよりは退化し ており、シミュレーションプログラムで高い性能を得るための手間が加速度的 に増大しているのが現状である。
特に話を面倒にしているのが、「分散メモリ型並列計算機の上でのメッセージ・ パッシングプログラミングモデル」というものである。これは、要するに、並 列計算機といってもそれぞれ独立なメモリと CPU をもったたくさんの計算機 の集まりであり、それらが互いに通信できるというだけのものである。これは もちろん並列計算機では非常に古くから使われてきたプログラミングモデル、 というよりは「計算機を作る側ではなにも考えないから、使うほうでなんか考 えてよ」という、計算機工学の側の仕事放棄を表しているだけともいえる。
もちろん、そのようなモデル、というよりはモデルの不在が特に科学技術計算 の分野で主流になったのは、使う側がそれを選択したからという面がある。問 題が単純でこのようなものでもなんとかプログラムが書けてしまうような場合 には、このようなやり方で高い性能が得られる。また、各プロセッサでは普通 の Fortran や C で書いたプログラムが動いているだけなので、プログラムを 書き始めるのに対する敷居はそんなに高くない。さらに、チューニングに関し ても、基本的には従来の単一プロセッサの場合と同じなので特別なツールとか がいるわけではない。これに対してもっと高度な並列処理言語のようなものを 使った場合には、多くの場合にプログラムの実行効率が下がるし、また高水準 のプログラムが実際にどのようにハードウェアにマッピングされているかが隠 されてしまうためにチューニングが容易ではなくなる。
とはいえ、ツリー法のように原理的にはすべてのステップが並列実行できる比 較的素性のいいアルゴリズムでも、これを分散メモリ型の並列計算機に載せよ うとするとアルゴリズムのかなり根本からの変更が必要になる。以下、そのあ たりを少し議論しよう。
ツリー法や FMM を分散メモリ型計算機に実装するうえで問題になるのは、空 間分割をどのようにやるかということである。ツリー法は基本的に遠くのほう は適当に計算するという方法なので、物理的に近くにある粒子を同じプロセッ サに割り当てることでプロセッサ間でのデータ交換を少なくすることができる。
ここで、粒子の分布が空間的にほぼ一様であれば、基本的に計算領域を規則的 に同じ大きさに分ければ粒子数は均等になるし、またその結果計算量もおおむ ね均等に分配されて高い実行効率が期待できる。が、広い領域に銀河が1つと か2つだけあって、まわりにはほとんどなにもないという状況ではそうはいか ない。
現在までに作られたほとんどのプログラムで使われているのが、再帰的直交分 割 orthogonal recursive bisection という方法である。 [WS92]この方法の原理は以下のようなものである。まず粒 子をある座標軸方向に並べて真中で分ける。分けたものそれぞれについて、今 分けたのとは直交する座標に並べてまた真中でわける、、、というのを、プロ セッサ数になるまで繰り返す。実際に並べなくても真中あたりで分ければいい わけなので、実装はいろいろである。
こうして出来た分割を各プロセッサに割り当て、基本的には各プロセッサで独 立にツリー構造を作る。
各プロセッサ毎にツリーが出来れば、自分のなかでの相互作用は普通と同じよ うに計算できる。別のプロセッサにあるツリーからの力については、まず、 「そこからの力を計算するのに必要になりえる範囲のツリー構造」を送っても らって、そこからの力を普通の方法と同じように計算すればいい。このツリー 構造は、基本的にはツリーの枝のあちこちを切り落としたものになる。つまり、 ある節点のところで、すでにツリーを受けとる側のどこからみても十分に遠 いという条件を満たしていれば、そこから下は決して計算に使われないので送 る必要がなくなるのである。
もっとも、この方法では若干の無駄がでる。というのは、本来ならば複数プロ セッサの粒子をまとめて1節点として重力が計算出来るような場合でも、そ のような大きな構造を作っていないのである節点からの力は少なくとも1度 は計算されることになるからである。また、 ORB の境界はツリー自体のノー ド境界と無関係なので、本来同じ節点に属する粒子が2つのプロセッサにま たがる結果計算量が増えるという効果もある。
これらがどれくらい大きな損失であるかに関するきちんとした定量的な評価は 残念ながら文献にはないが、プロセッサ数が非常に大きくなると問題になるの は明らかである。例えば、プロセッサ数が数百程度になると、プロセッサ毎 に一つというのでも本来のツリーで必要な計算量よりも大きくなる。また、通 信についても、データ量はともかく回数が増えるためのオーバーヘッドが無視 できなくなってくる。
この問題を回避する方法は2種類知られている。一つは通常の8分木をやめて、 各プロセッサの中で ORB をさらに繰り返し適用することでツリーを作るとい う方法である。8分木に比べてツリーを作るコストはかなり増えるが、まあ通 常の実装ではもともとそこの計算コストは小さいので大した問題ではない。
もう一つは、8分木を使って、なんらかの方法で全体ツリーを作り直すことで ある。原理的にはベストな方法は、各プロセッサが他の全プロセッサから必要 な情報をもらってきた後で、その情報を使ってそのプロセッサの粒子からみた 全体ツリーを作るというものである。このためには、単に他のプロセッサから 渡ってきたツリーの葉になっているものを全部粒子だと思って、もともと自分 が持っている粒子とまとめてツリーを作り直せばいいことになる。
まあ、これは少なくとも自分のが既に持っている分については2度手間になる。 普通は既にいったようにここの部分の計算時間は知れているので問題ではない が、ここの部分が問題になるようであれば減らすことは出来る。それには、ツ リーを作る方法を変えればいい。
基本的なツリー法のアルゴリズムを説明したところでは、ツリーを作る方法と して最初に全粒子が入る大きな立方体を考えて、それを分割していくという方 法を紹介した。これはまあわかりやすい方法であると著者は思うが、 Barnes と Hut の最初の論文では全く違う方法が説明されている。先に説明した方法 は、実は著者がベクトル化を可能にするために導入したもの[Mak90]である。
Barnes と Hut の方法では、まずからっぽの箱を作り、そこに一つずつ粒子 を入れていく。ある箱に入っている粒子が2つになったら、その箱を分割して 粒子を振り分ける。分割してもまだ同じ箱に入る場合ももちろんあり、その時 は単に分割を繰り返す。
このアルゴリズムが実装されていれば、他のプロセッサから渡ってきた情報を まとめるには、それらの「粒子」を既にあるツリーに入れていくだけである。
余談になるが、ツリー法を共有メモリ型の計算機で並列化するという研究はか なりたくさんある。これらの多くは、計算機科学（並列処理）の研究者によっ てなされたものであった。その多くは ツリーを作るのに Barnes と Hut のも ともとの方法を使って、それを非常な苦心をして並列化している。代表例は良 く知られた Stanford 大学の SPLASH ベンチマーク[WOT+95]である。 並列化するには複数の粒子を同時にツリーにいれたいわけだが、一つの粒子の ために構造を変えるのと他のどれかのために構造を変えるのが干渉すると結果 が変になってしまうので、同期とかロックとか難しいテクニックを使って矛盾 が起きないようにする。もちろん、その結果として速度は低下し、ある程度以 上プロセッサが増えても速度は上がらなくなる。これに対し、著者の方法は、 ベクトル化も共有メモリ型の計算機での並列化もトリビアルであり効率も良い。 なにがいいたいかというと、並列処理の研究者の場合、どうしてもアルゴリズ ムを所与のものとしてそれと「同じ」計算手順を並列に実行しようという方向 に目が向く。このために、結果が同じでもっと簡単に並列化でき、さらに並列 実行でない時でも性能差が問題にするほどはないような簡単な方法が作れる場 合でも、そういったものを作るということにそもそも考えがいかないことが多 いようである。
思いきり話がそれたが、このようにして全体のツリー構造を各プロセッサで作っ てしまえば、そこから先の各粒子への重力の計算は各プロセッサで全く独立に 行なうことができる。このために、ツリー法の並列化の効率はきわめて高く、 プロセッサ数が 1000を超えるような場合でも台数に比例した性能向上を実現 した例がいくつかある。
さて、今までは、ツリー法の場合、時間刻みは全粒子同じで一定である、つま り、独立時間刻みは使わないという前提で話をしてきた。銀河全体のシミュレー ションや、あるいは銀河団、宇宙論的構造形成の計算のように、一定の時間刻 みとポテンシャルソフトニングを組み合わせるのが効率的な応用がいろいろあ るのは確かである。しかし、球状星団のシミュレーションのようにそれではす まない応用もまたたくさんある。ここでは、そういった応用向けに、ツリー法 と独立時間刻みを組み合わせる方法について検討する。
といっても、並列化、特に分散メモリの計算機での並列化を考えなければ、こ れはそれほど難しい話ではない。基本的には、ある粒子への重力を計算するの に、他の粒子すべてからの力を合計するわけではなくツリーをたどって計算す るというだけなので、ツリーの節点に対しても予測子を持たせればいいわけ である。重心運動については子節点（あるいは粒子）の予測子から容易に作 れるし、4重極も式は繁雑だが計算するだけである。ある粒子を動かした時に は、その度にその粒子の上の節点の予測子はすべて更新する必要がある。 とはいっても、手を抜くことはできなくはなくて、実際に予測子が呼ばれるま で計算しなおさないというふうにすれば多少速くなる。
もう一つ重要なのは、木構造が変わる時にどうするかということである。粒子 が動き、自分のもともとはいっていたセルから出ていくので、それをほおって おくとセルの幾何学的大きさが粒子分布を反映しなくなる。このために、幾何 学的な判定条件を使って重力計算すると誤差が大きくなってしまう。 実際の粒子の広がりを計算しなおして、そっちを使って判定することで計算誤 差の増加は防げるが、もともとのセルサイズよりも粒子の広がりは大きくなっ ていくので計算量が増える。従って、適当なところでツリーを作り直す必要が ある。
もっとも簡単な方法は、適当な間隔で木構造をゼロから作り直すことである。 宇宙論的シミュレーションや銀河形成のための計算コードでツリー法と独立時 間刻みを組み合わせたものではほとんどこの方法がとられている。
これに対し、ツリー構造が変わったところだけ組み直すというのは実はそれほ ど難しくはない。これにもいくつかの方法があるが、一つは、ある粒子がある セルから別のセルに動いた時に、その2つのセルの共通の祖先であるセルの下 は全部作り直すという方法である。これは McMillan と Aarseth によって実 装された。[MA93a]もう一つは、実際にツリーから一旦問 題の粒子を取り除き、もういちど付け加えるという方法である。著者はこの方 法を実装したことがあって、たいして面倒ではないが現在のところ広く使われ てはいない。
現在のところ、ツリーと独立時間刻みとの組み合わせはシミュレーションのカ バーする時間範囲（系のダイナミカルな時間スケールで測った）が非常に短く、 そのかわり扱う粒子数が多い宇宙論や銀河形成のような問題、あるいはあまり 空間構造が発達しない惑星形成のシミュレーションのような場合に限って使わ れている。球状星団における自己相似的な密度上昇のような、系のわずかな部 分だけが極端に短い時間スケールを持つという場合に、どのように並列化をす ればよいかがまだほとんどわかっていないからである。
可能性としてはうまくいくかもしれない方法としては、 Springel らが導入し た並列化の手法がある。[SaSDW00]彼らの方法では、前に述べ たように ORB で空間分割して粒子をばらまくが、あるタイムステップで動か す粒子への力を計算するときに、その粒子があるプロセッサが他のプロセッサ から必要になる情報をもらってきて重力を計算するのではなく、逆に動かす粒 子の座標の方を他のプロセッサにおくって力を計算してもらい、力を送り返し てもらって合計するという方法をとる。力の計算が分散しているぶんだけ、粒 子があるプロセッサがその粒子への力の計算を全部やる方法に比べれば効率が 上がる可能性がある。
とはいえ、もともと独立時間刻みを効率的に並列化するだけでわりあい面倒な ので、ツリー法が付け加わるともっとうまくいかなくなるのはある意味やむを えないところがある。
結局のところ、並列化がうまく出来て効率をあげられるためには、もとの問題 に十分な並列性がある必要がある。ところが、独立時間刻みにしてもツリー法 にしても、単純なやりかたに比べて計算量を減らしてはいるがその分だけ並列 性も減っているわけである。
もっとも、粒子数があまり大きくない問題に対してツリーと独立時間刻みを組
み合わせる方法の研究があまり進んでいない理由はもうすこし違うところにあ
るかもしれない。一つは、現在扱える程度の粒子数 (数万以下)に対してはツ
リー法よりも前に述べたネイバースキームのほうが有利であり、またこちらの
ほうが並列化についても（おそらくは）容易であることである。粒子数が数万
程度の時、ネイバースキームでは計算精度をほとんど落さないで直接計算の
10-20 倍程度加速できる。これに対してツリー法では計算量 1/10 に程度にす
るには精度をかなり落す必要がある（具体的には力の精度で 
程度）。
それで得た答が正しいかどうかが良くわからないので、あまり使いたくないわ
けである。
もう一つは、専用計算機の存在である。これについては次章でもうすこし詳しく 述べる。
というわけで、ここまでで重力多体系の数値計算を普通の汎用計算機でやるという話は 一応おしまいである。
