というわけで、この講義も最終回である。予定では
Goodman, J., Heggie, D. C. and Hut, P., On the Exponential Instability of N-Body Systems , APJ 1994, 415, 715.
Fukushige, T., Makino, J., Nishimura, O. and Ebisuzaki, T., Smoothing of the Anisotropy of the Cosmic Background Radiation by Multiple Gravitational Scattering , PASJ, 1995, 47, 493-508.
Funato, Y., Makino, J., and Ebisuzaki, T., Attenuation of luminosity of distant galaxies by multiple gravitational lensing , APJL, 1994, 424, L17--20.
Chernoff, D. F. and Weinberg, M. D., Evolution of globular clusters in the Galaxy , APJ, 1990, 351, 121--156.
Fukushige, T., Heggie, D. C., Pre-collapse evolution of galactic globular clusters , MNRAS, 1995, 276, 206--218.
これだけでおわりにしては何なので、少し役に立つ（かもしれない）話をして おく。
ここまで、自己重力多体系での物理過程の話をしてきたわけだが、では、いっ たいそういった系をどうやって研究するかという話になる。球対称な系とか、 いろいろ仮定をおいていい場合には、無衝突系の場合には無衝突ボルツマン方 程式を直接数値積分するというような方法もあり得る。また、衝突系の場合に は orbit-averaged Fokker-Planck 方程式を数値計算するということも考えら れる。
しかし、次元を落した系でなければ、上のような方法は現実的ではない。計算 に必要なメモリ量や計算時間があまりに莫大なものになるからである。このた め、どちらの場合でも、実際に多粒子系の数値計算を直接行なうという方法が とられてきた。
解くべき問題は重力多体問題、つまり重力で相互作用する質点系がどのように ふるまうかという問題である。
系の進化を表す方程式は
つまり、単にある粒子はほかのすべての粒子からの重力を受けるというそれだ けである。
解くべき常微分方程式系は式10.1で与えられるので、初期条件を作 ればあとは各粒子の位置からそれぞれへの重力を計算するサブルーチンを書き、 それを適当な数値計算ライブラリに渡して数値積分すればそれでおしまいと行 けば話は簡単だが、それでは今日わざわざこの講義をする必要はない。そうう まく行かないのには2つの理由がある。
第一の理由は、ナイーブな方法では1ステップ当たりの計算量が粒子数の2乗に
比例して増えてしまうということである。もちろんもっとうまい方法がなけれ
ばしょうがないが、残念なことに計算量を 
から 
ない
しは 
に落とす方法が知られている。これらの原理は、基本的には遠く
の粒子からの力は適当にまとめるということである。
もう一つの理由は、すべての粒子を同じ時間刻みで積分すると、ステップ数が 非常に大きくなってしまうということである。例えばたくさんの微惑星が太陽 のまわりをまわっているという計算をしたとする。ほとんどの粒子に対しては、 太陽からの重力が圧倒的であり、また軌道は円軌道に近い。従って、例えば一 周期を100ステップとかそれくらいで刻んでやれば十分正確に軌道を追いかけ ることができる。しかし、いくつかの微惑星は他の微惑星と近接遭遇をしてい るかもしれない。するとこの近接遭遇をちゃんと計算するためには時間刻みを 短くする必要がある。このときに、すべての粒子の軌道を同じ時間刻みで積分 していると、別に近接遭遇中でもなんでもない粒子もすべて短い時間刻みで計 算されることになり、ステップ数が非常に大きくなる。この問題は、粒子毎に バラバラな時間刻みを与え、必要に応じてそれらを独立に変化させることがで きれば解決できる。
なお、第二の問題、すなわち近接遭遇の問題には、別の対応もある。それは、
相互作用のポテンシャルを変えて、近接遭遇が起きても時間刻みを短くしなく
てもいいようにすることである。具体的には、ふつうの重力ポテンシャルは
で与えられるが、これを例えば
といった
形にしてしまう。このようにすれば、 
より小さいところでは重力
が大きくならないので、粒子の典型的な速度を v とすれば時間刻みを
にくらべて十分小さく取っておけばいい。これをソフトニングと
いう。
ソフトニングをしていいかどうかは問題による。例えば微惑星の計算では、近 接遭遇で軌道が変化する、あるいは物理的に衝突して合体するといった過程を 計算したいので、ソフトニングをしたのではなんだかわからなくなってしまう。 これにたいして銀河のシミュレーションとかいった場合には、近接遭遇で軌道 が変化するといった効果は非常に小さく、ソフトニングが有効となる。
銀河、銀河団のシミュレーションでは、ソフトニングが使えるためにタイムス ケールの問題は無い。このために微惑星とか球状星団にくらべれば話は簡単で あり、第一の問題、すなわち計算量が粒子数の2乗に比例して増えるという問 題だけ考えればいい。
とはいえ、これは大きな問題である。例えば普通の渦巻銀河を、２体緩和によっ
て壊れてしまわないように表現するためには少なくとも
個程度の粒子が
必要である。このような銀河 1000 個から銀河団を作れば粒子数は 1億個とな
り、1ステップ計算するのに 
演算程度が必要になる。これは例えば
100 Gflops の超並列計算機で10日程度であり、 1万ステップ計算するのに300
年かかってしまう。
もう少し速く計算するためにひろく使われているのが Barnes-Hut アルゴリズ ムあるいはツリーコードと呼ばれる方法である。原理は簡単で、遠くのほうの 粒子はグループにまとめ、グループからの力は多重極展開（といっても普通に 使うのはせいぜい1次または2次、1次の展開は要するに重心で置き換えること） で計算する。
グループ分けは粒子毎に作っていてはしょうがないので、空間を階層的な立方 体からなるツリー構造に組織してあらかじめすべての階層で多重極展開を計算 しておく。ある粒子へのある立方体からの力は、それらが「十分に離れて」い れば多重極展開を計算して得られる。そうでなければ下の階層に再帰的に降り ていって、十分離れているとみなせるところまでいって計算する。系全体から の力も、単に系全体に対応する立方体からの力を計算することで求められる。
この方法では計算量が 
から 
になる。もう少し（原理的
には）賢い方法として、 
になるアルゴリズムが知られている。この方
法では、力を及ぼすほうだけでなく受ける方もまとめて計算する。つまり、力
を及ぼすほうの回りで有効なポテンシャルの多重極展開を、受ける方の中心で
のテイラー展開に変換し、そのテイラー展開を各粒子の位置で計算することで
力が求められる。
の計算法は天文シミュレーションでは広く使われている。また
パイプライン型ベクタプロセッサや並列計算機での実行もかなりよく研究され
ている。これに対し
のアルゴリズム (Fast Multipole Method, FMM)は
まだ天文の問題に使えるような実用的な実装（実際に走るプログラム）を作っ
た人がいないので、以下ツリーコードについてその実装上の問題点を簡単にま
とめる。
ツリーコードでどのような問題が発生するかを考えるためには、比較の対象と
して 
の直接計算アルゴリズムを考えることが有用である。これは、
ベクタプロセッサでは非常に高い効率で走る。またほとんどいかなる並列計算
機でも、マシンの理論ピークに近い性能を出すことができる。
もっともピーク性能を出すのが困難であると考えられる分散メモリ MIMD 型の 時にどのような計算法を使えるかを簡単に説明すると以下のようになる。
まず各プロセッサに（大体）同じ数の粒子を割り当てる。自分ので持っている 粒子同士の相互作用はそのまま計算すればいい。それ以外の他のところとの総 誤差用は、プロセッサを1次元的なリングとみなして、粒子をたらい回しにし て順次計算することができる。また、あるプロセッサのデータを全プロセッサ に効率良く放送できる場合にはその機能を使ったほうがプログラムが容易であ る。
計算速度に対し必要な通信速度は、プロセッサ数をpとしてオーダとしては
になるので、p<<N（普通成り立つ）では通信速度は問題にならない。
共有メモリならば話はもっと簡単になる。
さて、ツリーコードの場合、演算量は減るが計算速度（1秒に何回浮動小数点 演算をするかという意味での）は、普通のスカラープロセッサでも低下する。 これは、ある粒子への力を計算する時の配列のアクセスの仕方が不規則であり、 キャッシュミスの確率が高いことが大きな理由である。また、そもそも分岐自 体が多いことも性能低化の理由となる。
ベクタプロセッサではより困難な問題、すなわち再帰呼び出しがあるようなルー プはベクトル処理できないという問題が発生する。これは87年頃にいくつかの アプローチによって解決された。 （詳細は省く）
これらのアプローチは、アルゴリズムの変更による計算量の増大またはメモリ アクセスの不連続性による計算速度の低下というペナルティを持つものの、もっ とましな方法は知られていないということもあり現在も広く使われている。
共有メモリ MIMD 型の並列計算機の場合、再帰呼び出しも並列化可能で特に問 題なく実行可能となる。ただし、分散共有メモリ（仮想共有メモリ）の場合、 物理的には分散メモリなので後に述べる分散メモリのマシンと同じ配慮が必要 になる。
分散メモリでツリーコードを実行する方法については、Caltech のグループは 87年頃からいくつかの方法を試み、まあまあの結果を得ている。基本的なアイ ディアは、粒子を空間的に近いグループに分けるということである。分けるの には、それ自体にツリーによる空間分割を利用する方法や、適当に空間を半分 に切っていく方法などが知られている。
第二の問題、すなわちタイムスケールの問題には、粒子毎に別々の時間刻みを 与えることで対応することになる。この場合、時間積分は 「次の時刻」すなわち現在の時刻と時間刻みの和がもっとも短い粒子を選び、 それの軌道を積分してその時刻と時間刻みを改めて計算し、最初に戻って次の 粒子を選ぶという具合にすすむことになる。
この方法の場合、ベクタプロセッサでの実行は容易である。というのは、ベク トルの総和というのはすべての現存のベクタプロセッサで効率良く実行できる 処理であるからである。
並列計算機の場合、一旦各プロセッサで自分の分を計算し、その総和をとって から軌道積分に進むということになる。これは p<<N であれば不可能ではな いが、時間刻み共通の場合に比べて難しい。これは、一回の通信でのデータ転 送量が小さいため、ピーク転送速度ではなくメッセージレイテンシ（最初のデー タを送り出してからつくまでの時間）が問題になるからである。実際に、必ず しも高い性能が得られないことがわかっている。
「独立時間刻み」法の改良として、「独立2重時間刻み」あるいは
neighbor scheme といわれる方法がある。これは、ある粒子への力をその回り
の粒子からの分と遠くからの分に分割し、遠くからの分は長い時間刻みで積分
するという方法である。理論的には、単純な方法に比べた時の計算量の減少が
に比例するということになっている。
この方法の並列化は困難である。近くの粒子からの力の計算ではメモリアク セスが不規則でしかもループ長が短くなるので、ベクタプロセッサでも効率が 大きく低下することが知られている。
並列計算機上で効率をあげるには、一度に1つの粒子しか動かさないのではな
く、なるべくたくさんの粒子を並列に積分できることが望ましい。これは時間
刻みを2のベキ乗に離散化することで可能になる。このようにすれば、ある粒
子が積分される時にそれと同じかあるいは短い時間刻みを持つ粒子は並列に積
分して良く、並列度は劇的に向上できる（理論的には 
と見積もら
れている）。現在までの分散メモリ型並列計算機への独立時間刻み法の実装は、
すべてこの方法（ブロック時間刻み）によっている。
ここまでは、「相互作用をいつ、どうやって計算するか」という話をしてきた。 しかし、実際に数値解を求めるには求まった重力加速度を使って各粒子の軌道 を積分する必要がある。以下、その方法について簡単に述べる。
これは知らない人はいないと思いたいが、ある微分方程式の初期値問題
があるときに、時刻 
での近似解 
を、
時刻 t での近似解
を使って
とする方法である。局所打ち切り誤差は 
であり、したがっ
て大域誤差は 
ということになる。プログラミングは簡単である
が、誤差があまりに大きくほとんど実用にはならない。
これはセパラブルなハミルトン系、すなわちハミルトニアンが運動量の関数と 位置の関数に分かれるものに対する特別な方法である。通常の、運動量に対す るハミルトニアンが2次形式である場合には、
と書くのが普通である。これでは速度と位置がずれた時間でしか定義されない が、出発用公式として
を使い、さらに終了用公式として
を使うことで最初と最後を合わせることが出来る。この形は、実は
と数学的に等価である（証明してみること）また、以下のようにも書ける
さらにまた、速度を消去して
, 
, 
の関係式の形で
かいてあることもあるかもしれない。なお、すべて違う名前がついていたりす
るが、結果は丸め誤差を別にすれば完全に同じである。
この公式は局所誤差が 
、大域誤差が
である。
さて、局所誤差という観点からはこれは決して良い公式というわけではないが、 現実には特に無衝突系の計算ではこの公式は非常に広く使われている。衝突系 の計算で、独立時間刻みを使うような場合にはこの方法は使われないが、それ 以外の非常に多くの問題はこの方法で解かれているといってよい。
これは、別にもっとよい方法を知らないからとかではなく、実はこの方法がい くつかの意味で非常に良い方法であるからである。いくつかの意味とは、例え ば力学平衡の系でエネルギーや角運動量が非常によく保存するということであ る。これらの保存量については多くの場合に誤差がある程度以上増えない。
この、ある程度以上誤差が増えないというのは目覚しい性質である。通常の方 法では、エネルギーの誤差は時間に比例して増えていく。従って、長時間計算 をしようとすればそれだけ正確な計算をする必要がある。ところが、エネルギー の誤差は溜っていかないのならば、かならずしも精度を上げる必要はないとも 考えられる。
もちろん、エネルギーが保存していればそれだけで計算が正しいということに はならない。が、理論的には、いくつかの重要な結果が得られている。まず、
Symplectic method については、以下のようなことが知られている
というわけで、 leap frog が、それよりも局所的には誤差が小さいルンゲ・ クッタとかに比べて良い結果を与えるのは当然のことである。
さて、独立時間刻みを使った場合には、話はそう簡単にはいかない。というの は、 symplectic method の良い点も symmetric method の良い点も失われる からである。このため、通常はより高次の方法が用いられる。まあ、 cosmological な計算なんかでは独立時間刻みで2次の公式を使っている例もな いわけではないが、結果を信用できるかどうかには注意したほうがよいかもし れない。
高次の方法としては、伝統的には Aarseth が60年代に開発した方法が用いら れてきた。これは、いわゆる線形多段法 linear multistep method に基づく ものである。具体的な実現についてはここでは省くが、基本的に、過去数ステッ プの加速度から補外多項式を作り、それを積分して位置と速度を更新する。独 立時間刻みを使うに当たっては、他の粒子の位置は補外多項式で計算する。
Aarseth method は広く使われてきたが、結構面倒であるし高次にしていくと 計算量が増え、また安定性が急激に悪くなる。もう少し簡単な方法はないかと いうことで考えられたのが Hermite method である。これでは、加速度の他に その時間導関数まで直接計算し、それから補外、補間多項式を作る。これは最 近かなり広く使われるようになってきた。
以下のどれでもOK。
