next up previous
Next: 4 テスト計算 Up: The Pseudoparticle Multipole Method Previous: 2 ツリー法の原理

3 PM法

高次の項を扱うのに GRAPE を使う方法はないかと考えると、「多重極展開を粒子で表 現する」ということができればよい。もともと多重極展開は実際の粒子分布の 作る場を近似しているわけだから、適当に粒子(以下、こちらを本当の粒子と 区別するために疑似粒子と呼ぶ)をおくことで多重極展開を表現 しなおすということは原理的には可能なはずである。この時に、元の粒子より も少ない数での疑似粒子で、ある次数までは元の粒子と同じ多重極展開を持たせる ことが出来れば、あとはそれらの疑似粒子からの重力を普通の粒子と同じよう に計算することでツリー法が実現できることになる。

重力の場合、いくつかの粒子の重心は p=1 まで展開が合うという意味でそ のような粒子のおき方の例になっている。

次に4重極項を考えてみる。重心回りの4重極項モーメントテンソルはトレース レスなので、対角化すれば独立な要素は2個である。言い換えれば、疑似粒子 は第一軸と第二軸の張る平面上におけばよく、3個の等質量の粒子を2等辺3角 形におくことでで4重極項は表現出来る。

もっと高次の場合でも、このように粒子配置を構成するのも原理的には可能なはずだが、 しかし、次数をあげていくと解くべき方程式がどんどんややこしくなるし、だ んだんどうやって式を立てればいいかも良くわからなくなる(わかる人は是非 やってみて欲しい)。従って、以下ではもうちょっと安直なアプローチを考え る。

結局やりたいことはなんであったかというと、多重極展開の係数が実粒子の分 布と等しくなるような疑似粒子の分布を見つけだすということである。ここで、 疑似粒子をおく位置を球面上に限ってしまうと、疑似粒子の方の分布の 多重極展開は、もちろん疑似粒子の質量分布の球面調和関数展開そのものであ るから、別にこういう制限をつけても疑似粒子のおき方を決めることは依然と して可能である。もちろん、位置の自由度のうちの1つを殺しているので、そ の分たくさん粒子がいることになるが、とりあえず大した問題ではないと信じ よう。この時には、結局球面調和関数展開の逆変換で球面上の質量 分布が求まり、それを適当に離散的な質点で近似することで疑似粒子を決める ということになる。

さて、依然として問題であるのは、どうやって球面上に粒子をおくかである。 これは、上で求まった球面上の質量分布の、なんらかの意味での最良近似であ るようなものを有限個の粒子で表現するという問題になる。ここで注意して欲 しいのは、粒子一つの多重極展開は位置に対しては非線形な関数であるが、質 量に対しては線形な関数であることである。従って、位置を固定して質量だけ を変えるようにして、 個の粒子を適当において質量だけを変えてやるこ とで、必ず 個の展開係数の値を合わせることが出来るはずである。

が、落ち着いて考えてみるとこれはまだいろいろ面倒そうである。計算の順番 が、まず展開係数を求めてから、その展開係数を満たすような質量分布を決め るということになる。もちろん後半の計算では解くべき線形方程式の係数はい つも同じなので、LU分解しておくか逆行列を求めておいて、それから計算すれ ば一応はいいことになる。ここでは、点の数が増えることはあまり気にしない で定式化やプログラムが簡単になる方向で考える。

形式的には、するべきことは実粒子分布を多重極展開し、その係数を逆球面調和関 数展開して仮想粒子の質量に戻すということである。仮想粒子の数が無限に多 い連続極限では、形式的には逆変換の式は極めて簡単で

と書ける。 ここで は球面調和関数展開の係数であり、関数 は球面調和関数である。 簡単になるのは球面上で球面調 和関数が完全正規直交系を作っているからである。有限個の仮想粒子ですます 時も、もしもその上で球面調和関数を有限次数で切ったものが直交系 をなしていれば、上の式と本質的に同じ簡単な変換式がつかえることになる。

p次までの球面調和関数が有限個の仮想粒子の質量が作る空間の上で 直交系をなすためには、その空間のなかで内積演算が定義されてそれが直交性 を満たしていればよい。そのためには、結局 2p 次の多項式の積分が厳 密に求められればよい。そのように仮想粒子を置く方法はいろいろありえて、 例えば 方向は p次ルジャンドル多項式の零点に置き、 方 向は 2p+1 個を等間隔に置くというものでも実現できる。が、この場合、仮 想粒子の重みづけが の値によって変わる。

もうすこしましな置き方はないかと思うわけだが、これは実は非常に昔から計 算幾何学のほうで研究があり、 spherical designs という名前で知られてい る。これはまさに我々が必要としている、

という性質を満たすものである。つまり、 t次の spherical design という のは数学的には以下のように定義される。

単位球面S上の座標点 の集合で、t次以下の任意の多 項式に対して

が成り立つようなもの。

これには古典的な研究[10]と、最近になって詳細に評価しなおしたもの[6]があって、 若干結果が違うが20次程度のものまで求められている。また、座標値そのもの は WWW 上gifに公開されて いるので、ファイルに落すのは非常に簡単 である。次数 p が小さい時には、正多面体の頂点といったものになってい るが、それ以上では基本的には切頭立方体の組合せといった方法で構成されて いる。点の数は結局ルジャンドル多項式を使う場合とあまり変わら ないが、数値積分の重みがすべての点で等しいためにプログラムが簡単で数値 的にも性質がよい。

さて、 spherical design を使うことに決めると、結局実粒子から仮想粒子へ の変換は以下のように計算できることになる。

まず、球面調和関数展開の係数 は、

で与えられる。ここで は粒子 i の質量であり、 はその位置の極座標表示である。関数 は球面調和関数であり、

で与えられる。ここで はルジャンドル陪関数である。 これと同じ展開係数をもつ半径 a の球面上での質量分布

を満たす必要がある。球面調和関数の直交性から、この質量分布

で与えられることになる。ここで積分を先ほどのK 個の仮想粒子上での数値 積分で置き換えれば、仮想粒子 j の質量が

であたえられることになる。ここで、展開が p 次で打ち切られていること に注意してほしい。

さて、これで必要な変換は与えられたが、いちいち係数を求めてからそれをも う一回変換するのは美しくない。実粒子の質量から直接仮想粒子の質量を求め ることを考えてみると、これは i, l, m の3重の総和になって一見計算 が大変なような気がする。しかし、実は三角関数と同様に球面調和関数にも加 法定理

というものがあり [ここで は方向の間の角度]、これを使えばm についての和は消えて、 結局

 

ということになる。

ここで注意してほしいことは、式(13)は実粒子から仮想粒子へ の変換だけでなく、ツリーの上のレベルに上がっていく仮想粒子から仮想粒子 への変換にも同様に使えることである。つまり、式(13)で必要 な変換はすべて表現されていることである。もうひとつ重要なことは、l の 総和の式は結局 l 次多項式なので、まとも に計算すれば 、漸化式を使えば , 区分近似多項式を準備して おけば で計算できるということである。

なお、ツリー法の場合にはこれで話は終りであり、あとはセルから各粒子への 力を計算するときに疑似粒子からの力の和として計算するというだけである。 従って、従来のツリー法にわずかな手直しをするだけで任意次数のものが 作れる。

FMM[4,5]についても同様なことが考えられるが、ここで は省略する。その理由は、現在までに発表されたほとんどの論文で、計算精度が 同程度ならばツリー法のほうが FMM よりも有意に速いという結果が示されて いるためである[3]。

ツリーは計算量が で FMM は なのだから、 FMM のほう が速そうなものであるが、実際にはそうはならない。ここではこの理由に深入 りはしないが、簡単にいうとツリーでは多重極を評価するところだけで誤差が入るのに対し、 FMM では多重極変換を局所変換に変換して、さらにそれ移動してというふうに 段階が増え、それぞれの段階で誤差が入ってくるためである。



Jun Makino
Wed Apr 14 20:49:22 JST 1999