重力多体問題は、原理的には単に大規模な常微分方程式の初期値問題である。 というわけで、普通に考えると、数値計算ライブラリパッケージに入っている ルンゲ・クッタ法や線形多段階法のサブルーチンをそのまま使えばいいような 気もするが、これで話が済むことはほとんどない。これは、既に述べた2つの 問題に対応しないといけないためである。
実際に、重力多体問題を扱う上で鍵となる手法は2つある。一 つは独立時間刻みであり、もう一つは安定な連星系をサブシステムと認 識し、それに対しては特別な扱いをする方法である。
独立時間刻みというのは、基本的には粒子毎にばらばらの時刻 
と
時間ステップ 
を与える方法である。時間積分は、以下の
ような順序で進むことになる。
が最小の粒子を選ぶ。
と、こう書くと、もちろん謎が発生する。軌道を積分するといっても、 それにはなんらかの方法で他の粒子からの重力を計算する必要がある。 つまり、その粒子の新しい時刻での他の粒子の位置が必要になる。しか も、例えばルンゲクッタのような方法では、現在の時刻、新しい時刻の 他に、その中間のいくつかの時刻でも他の粒子の位置が必要になる。
例えば、時間積分が前進オイラー公式のような1次精度のものであれば、 これをどのように解決すればいいかというのは自明である。つまり、別 に何もしなくていい。1次精度の範囲内では、ある粒子の新しい時刻での 加速度を計算するのに、他の粒子についてはそれらのそれぞれの時刻で の位置をそのまま使ってかまわないからである。
ここで、とりあえず満たすべきは積分公式の次数、すなわちタ イムステップの何乗で誤差が小さくなるかということだけであるとする と、とりあえず、以下のような方針が考えられる。
といってもコトバを知ってる人にしかわからないような気がするので、以下説 明。
まず線形多段階法とは、 常微分方程式の初期値問題の解法の一つである。ルンゲ・クッタとの違いは、 何ステップか前での導関数の値や数値解そのものをとっておいて、それらの線 形結合として新しいステップでの解を求める方法である。
原理的には非常にいろいろな方法があり得るが、現在実際に使われているのは
ほとんどアダムス型公式といわれるものでる。この方法では、新しい時刻 
(ここでは
添字は粒子の番号ではなく時刻を表すことにする)での解 
を、前の
数ステップでの導関数から作った補間多項式 
を 
から
まで積分して求める。つまり
補間多項式の作り方にも無限にいろいろありえるが、実際に使われているのは
と p+1点とって p 次多項式を作る方法
と、
の p+1 点をとって同じくp 次
多項式を作る方法である。前者の方法では、 
を作るのにそこでの
導関数の値は必要ないので、これは陽解法ということになる。後者では 
を求めるのにそこでの導関数が必要になるので、これは陰解法である。陽公
式の利点はプログラムが簡単でステップ当たりの計算量も少ないことであり、
陰公式の利点は安定性、誤差の両方がステップサイズが同じなら
陽公式に比べて圧倒的によいということである。陽公式はアダムス・バシュ
フォース公式、陰公式はアダムス・モールトン公式という名前がついている。
余談になるが、アダムスは19世紀末のイギリスの数学者（物理学者・天文学者） であり、彗星の軌道の数値計算のためにこの方法を開発したそうである。
実際には、陽公式と陰公式を組み合わせて使う方法が実装も容易で便利である。つまり、
最初に陽公式を使って 
を推定し、それを使って導関数を計算し、
計算した導関数を使って今度は陰公式を適用して
を計算しなおす。
原理的には、さらに導関数を計算し、それを使って
を計算しなおす
のを、
が収束するまで繰り返さないと陰的解法にはならない。が、
通常の場合には、反復を全くやらないでおしまいとか、あるいは一度だけやっ
ておしまいにする場合が多い。
これは単に用語の説明だが、反復を全くやらない方法が PEC モードというも
のである。これは、陽公式の適用のことを predictor step （P, 予測）、導関数
の計算を evaluate step (E, 評価)、陰公式の適用を corrector step (C, 修
正)と呼ぶことにして、それらを一度ずつこの順番でやるという意味である。
これに対して、1度反復を繰り返す方法は 
モードと書かれ
る。原理的には PE とか PECE 
とかいった、Eで終る方法も
あり得る。
普通の数値計算の教科書等を見ると、
や 
のような公式が勧められていることが多いが、これは（後に述べ
る対称型公式のような）特別な場合を除けば重力多体問題ではあまり得にはな
らない。これらのほうが、真の収束値に近い分ある程度誤差が減ることは確か
である。しかし、計算量のほとんど全てが加速度の計算にかかるので、タイム
ステップが同じなら計算量がほぼ2倍になる。それよりも PEC モードでタイムス
テップを半分にしたほうが得になることが多い。
なお、次数についてだが、これはもちろん粒子数、要求する精度等、さまざま な条件によって変わってくる。例えば太陽系の惑星の長時間積分というような 話になると、12次とか14次といった超高次の公式が使われる。これに対し、球 状星団や銀河中心核といった問題では4次程度の比較的次数の低い公式が使わ れることが多い。理論的にも実験的に要求する計算精度に依存して最適な次数 が存在するわけだが、現状では最適な次数というよりは実装の単純さ等によっ て次数を決めている面もある。
なお、この方法でプログラムを実際に書くのは結構面倒である。この後で、も
うちょっとプログラムが書けるような方法について議論するが、基本的な考え
方というか、キーワードは示しておこう。タイムステップと次数が固定の場合
は話は簡単である。アダムス法の場合は、結局 
が
導関数の値の線形結合で書けるわけなので、その係数をあらかじめ求めるか、
適当な数値解析の教科書から写してくればいい。
可変刻みにすると話が急にややこしくなる。線形結合で書けるといっても、そ
の係数が毎回変わるからである。実際のプログラムでは、まず不等間隔補間多
項式を構成し、それを代数的に積分（もちろん多項式なので各項を次数+1 で
割るだけだが）し、求まった積分された多項式にタイムステップをいれて
を求めるという手順になる。
補間多項式を構成する標準的な方法は2つあって、一つは差分商 (divided difference) を使ってニュートン補間多項式を作る方法である。これは、もち ろん求まる多項式は一意に決まるのでラグランジュ補間と同じだが、数値的に 安定であるとか、また補間範囲を一点ずつずらしていくときに前の計算結果の 使い回しによって計算量を減らせるとかいった利点があり、名前からも 想像できるように古くから使われている補間多項式の構成法である。
もう一つは、通常 multivalue method または提案者の名前をとって Nordsiek method といわれている方法である。この方法では、差分商の代わりに補間多 項式の係数、つまり高階の導関数の値そのものを記憶しておく。これの良いと ころは、見かけ上1段階法になっているので気分的にタイムステップを変える のが簡単であるということである。なお、実際には、見かけ上1段階法である が、補間多項式をずらすときには、考え方としては等間隔補間の元の参照点で の値を作り直し、もっとも古い点を落してそのかわりに新しい点を付け加えて 補間多項式を作りなおすというような手間が必要になる。タイム ステップが頻繁に変わる場合には、 multivalue method では補間多項式の性 質が多少悪くなるという問題が少なくとも理論的には存在する。
以下、もう少し積分公式の話を続ける。
まず、細かい話であるが、ハミルトン系のような2階の方程式の場合には、ア ダムス型公式をそのまま使うのはあまり賢い方法とはいえない。通常、高階の 微分方程式の数値積分をするには、最高次の一つ下までの導関数を変数にとっ て連立の一階微分方程式系に置き換える。が、アダムス法の場合、加速度に対 する補間多項式を2回積分すれば位置の変化が求まるわけであり、速度につい ての補間多項式を構成したり、速度の昔の情報をとっておく必要はない。この ため、必要なメモリ、計算量の両方が減る。
もう一つ、線形多段階法、あるいは一般に常微分方程式の数値解法のほとんど 全てについていえることだが、良く研究されている解法では基本的には微分方程式を 与えられたものと考え、その右辺値、つまり導関数だけを使って積分公式を構 成する。
多くの教科書で、高次の導関数を使って直接テイラー展開を計算して近似解を 構成する方法が「テイラー法」として紹介はされるが、通常高次の導関数を計 算するのに必要なコストは次数の指数関数で上昇するのであまり現実的な方法 ではないとしてそれ以上詳しい議論はなされない。
しかし、一般に非常に高い次数まで計算すれば指数関数的に高価になるのは確 かだが、例えば次の導関数くらいならそれほど大変ではないことが多い。例え ば、2体間の重力
なら、その一階導関数は
である。 ここで
であり、 
、
はそれぞ
れ粒子 i の位置、速度、Gは重力定数、
は粒子jの質量である。
式はちょっと複雑だが、追加計算には除算や平方根演算は必要ないの
で、計算時間の増加はたいしたことはない（計算機によるが、2倍以上になる
ことは稀である）。これ以上高階の導関数を計算するのは、さすがに現実的で
はない。
このように微分方程式の右辺の時間導関数が比較的安価に計算できるので、ア ダムス法を拡張して、補間多項式を導関数も使って構成することが考えられる。 この型の補間自体はエルミート補間と呼ばれる。エルミート補間を使う積分公 式は、人によってエルミート型、あるいはオブレヒホフ型の公式と呼ばれる。 この公式自体は古くから知られていて、数値解析の教科書を見ると載っている ものもある。が、まあ、一応、重力多体問題にこの型の公式を導入したのは著 者であるらしい。[Mak91b]
この方法の 利点は、段数に対して達成出来る多項式の次数が2 倍になることである。特に、 実質的には1段法である2段の陰的公式で4次精度が実現出来るため、多くの応 用ではこれで十分な精度が達成出来る。
ただし、この場合予測子は2段とはいえやはり多段階法であり、実装はけっこ う面倒である。ここで、さらにもう少し手を抜くために、予測子の最高次の次 数を位置の時間に対する3次までで止めるという方法が考えられる。誤差の係 数は別にして、公式の次数に関する限りは、これでもつじつまはあっている。
このやり方の大きな利点は、出発公式が不要であるということである。線形多 段階法一般にいえることだか、公式自体のプログラムよりも、最初に時間積分 を始めるところをどうするかのほうが面倒である。というのは、公式は加速度 の過去の値を必要とするのに、積分を始める時にはもちろんまだ過去の情報は 持っていないからである。
一般には、この場合にどうするかというと2通りの方法がある。もっとも一般 的な解決法は、タイムステップだけでなく次数も可変な公式を使うことである。 最初は次数2の公式から出発して、（必要に応じて短いタイムステップにして おき）、積分を進めて過去の値が入るのに応じて次数を上げていく。
次数が4程度と低い場合には、もっと直接的な方法もある。例えば、広く使われて いるプログラムでは、加速度の3次までの導関数を直接計算し、それから過去 の差分商を構成する方法が使われている。
というわけで、一般の線形多段階法で 2 次以上になると出発公式をどうする かを考えないといけないが、エルミート公式では 4 次でも出発公式が不要に なる。これは手間を省くという点からは馬鹿にならない。
以下、この場合に積分公式がどうなるかをまとめておく。
予測子は単に以下のようなテイラー展開で与えられる。
ここで 
と 
はそれぞれ位置と速度の予測子であり、
, 
, 
および 
は前の時刻での
位置、速度、加速度と加速度の一階時間導関数である。
は時間刻
み 
である。
この位置と速度の予測子を使って新しい時刻 
で加速度とその導関数
および 
を計算する。これらを使って、時刻 
での
加速度の2階と3階導関数は以下のように求められる。
これから、修正子を書き下して整理すると、
ということになる。なお、ここでは、最高次の項は速度の修正子にだけ使って、位置 と速度を同じ形に書き下した。最高次の項を位置に入れても、時間積分の精度 は上がらない。独立時間刻みにするなら、他の粒子の予測子も計算しない といけないが、これらも単にテイラー展開を評価するだけなのでなんというこ とはない。
あと、決めないといけないのは、タイムステップの計算式である。伝統的には、 以下のなんだか意味が良くわからない形の式が使われている。
これは、どういう理屈があって決まるのかといわれても困るような代物だが、 経験的には非常にうまく動くことが知られている。以下、「うまく」動くとは どういうことかをもうちょっと考える。
理屈としては、時間ステップをどのように決めるべきかというのは結構良くわ からない問題である。もちろん、位置なり速度の誤差の推定は出来るわけで、 その誤差をある値以下になるように調整しようというのが普通の話だが、ここ で問題になるのが空間構造が発達し、タイムスケールに大きな幅が出来るとい うことである。このために、例えば位置や速度の絶対誤差に上限をつけるとい うのがよい方法かどうかはあまり自明ではない。といっても、ではといって相 対誤差にすればよいというものでもない。つまり、絶対誤差にすれば、特に速 度が小さい粒子に対して相対誤差が非常に大きいということになる。ではといっ て相対誤差にすると、まず誤差がガリレイ変換に対して不変でないという問題 が生じる。例えば、系が全体として運動していると、そのなかでの粒子の運動 がどちら向きかで時間刻みが変わることになってこれは非常に変である。また、 もちろん、速度が大きいと絶対誤差としては大きくなるわけで、例えばエネル ギー保存を重視するなら、誤差のほとんどが速度の大きいところから出るとい うことになってなんとなく無駄なような気がする。
つまり、原理的には、多体計算において「どのような誤差を最小にするべきか」 という問題があるわけだが、これは、要するに、良くわかっていないのである。 もちろん、良くわかっている問題もある。例えば、太陽系の惑星の運動の長時 間積分という場合には、最終的に問題になるのは例えば 10 億年積分したあと でのそれぞれの惑星の位置そのものであり、それらの誤差が許容範囲内になる ようにタイムステップを決めることになる。また、逆に、非常に大きなスケー ルでの宇宙の構造形成といった場合にも、観測と比較するのは例えば銀河団の 空間分布の2体相関関数であり、さらに、その構造は自由落下のタイムスケー ルで出来たものである、つまり、いい方を変えれば、そのスケールでの構造は 一様な宇宙膨張から密度揺らぎが自己重力で成長することで出来るわけだが、 これはまだ成長中であって共動座標でみてものはあまり動いていないという事 情があるので、ここでも位置、速度そのものを見るというのはそれほど悪いや りかたではない。
問題は、この2つの中間の場合である。この領域では、シミュレーションの最
終状態でのそれぞれの粒子の位置は原理的にほとんど意味を持たない。これは、
各粒子の軌道そのものは、リヤプノフ指数が正で、その時間スケールが自由落
下時間 
(式85)の数分の一程度であるという意味でかなり
強くカオス的であるからである。これは、大雑把にいえば有効数字が
あたり1桁程度減っていくということであり、従って 10 
も
たてば初期条件の精度が 64 ビット浮動小数点数の有効数字だけあったとして
も、最終状態での個々の粒子の位置の精度は10桁落ちて 5-6桁しか残っていな
い。つまり、2つ全く同一な系で、1つの粒子だけ初期条件の最後の1ビットだ
け変えて計算を始めたとして、10 
程度でそのずれが
倍にも増幅されるわけである。
これは、力学系自体の性質なので、どれほど数値積分を正確に行なっても同じ ことである。球状星団や散開星団の進化の計算では、数千から数万クロッシン グタイム程度の長さのシミュレーションを行なうが、その後で各粒子の位置が 意味を持つためには計算精度が数万桁いるわけである。そんなことはできるは ずがないので、多体シミュレーションでは必ず個々の粒子の位置は「正しく」 ない。
このように指数関数的な成長が起きるのは、難しげにいうと
空間
（6N 次元位相空間）での軌道が双曲的であるからということになるが、も
ちろんこれは指数関数的に広がるというのの言い替えでしかない。具体的にミ
クロスコピックにみてどうやって広がるかというのを見てみると、これは結局
他の粒子と近接遭遇した時に軌道間の距離が拡大される傾向があるからである。
2次元の時のイメージを図6
に示す。
ここで、細かい人は 3次元だと紙に垂直な方向では軌道が縮んでいるの ではないかと思うかもしれない。軌道のずれの1次（線形）の分だけをとって 計算するとそうなるが、ここでは実は2次の項の分が効き、そちらは拡大する 方向に働くのである。
なお、ここで時々誤解[GS86]があるような気がするので
断っておくが、このように初期条件の差が指数関数的に成長する時間スケール
（大雑把にいって最大リアプノフ指数の逆数程度）と、系が熱力学的に緩和す
る時間スケールとの間にはなんの関係もない。実際、前に述べたように熱力学
的に緩和する時間スケールはほぼ粒子数に比例するのに対し、指数関数的成長
の時間スケールは粒子数に無関係に一定なのである。さらに、指数関数的成長
を担うのは、上に述べたように比較的近接した粒子同士の相互作用である。具
体的には、距離が 
の程度のものの寄与が最大である。この程度
の距離まで粒子同士が近付くのはクロッシングタイムに1度程度であり、これ
は粒子数 N に無関係である。これに対し熱力学的緩和を担うのは 
か
ら R までの全ての粒子との相互作用であり、従って近距離でのポテンシャ
ルの形をすこし鈍らせることで、熱力学的緩和の時間スケールをあまり変えな
いで指数関数的成長の時間スケールのほうは非常に長くすることができる。
ここでよろしくないのは、熱力学的緩和とは「系が初期条件を忘れるこ と」であり、また、初期条件の微小なずれが指数関数的に成長することがまさ にその「忘れる」ことであるとなんとなく思ってしまうことである。
少し思考実験をしてみよう。今、非常に粒子数が大きい系を考え、そのなかで、
ある位相空間の微小体積を考える。熱力学的緩和によって初期条件を忘れると
は、基本的にはその中にある粒子が系分布関数全体（というか、まあ、とにか
くある程度の広さ）に広がるということであろう。それらが依然としてどこか
に固まっていれば、緩和したとはいい難い。では、指数関数的な成長でこの微
小体積をどのように変形できるのであろうか？配位空間での大きさが、上の
よりも小さいうちは、空間の少なくとも1方向に向かってはほぼ
指数関数的に伸びていく。これは、他の粒子との散乱が基本的にはコヒーレン
トであるからである。つまり、微小体積の大きさに比べて散乱のインパクトパ
ラメータが大きく、微小体積のなかの各粒子の軌道は基本的には同じように曲
がる。同じように曲がるので、その差を考えると散乱するほうの粒子に近いほ
うが大きく、遠いほうは小さく曲がり、差が指数関数的に成長できるのである。
逆に、どんどん延びていってある1方向に対して長さが
を超える
と、散乱のインパクトパラメータに比べて微小体積のほうが大きいような散乱
が頻繁に起きることになる。このような散乱では、上の場合のように線形に微
小体積が引き延ばされるというような簡単な話にはならない。むしろ、長さが
よりも十分大きい極限では、各粒子がどのような散乱を受けるか
というのは他の粒子とはほとんど無関係になる。もちろん、非常にインパクト
パラメータの大きな散乱に対しては依然としてコヒーレントなわけだが、それ
による指数関数的成長の時間尺度はインパクトパラメータが大きくなるにつれ
て急速に長くなる。こうなると、各粒子間の距離は指数関数的には広がり得な
い。散乱が無相関である極限では粒子間の距離は拡散的に広がることになる。
そういうわけで、指数関数的に軌道が広がるといっても、それはせいぜい
程度の距離までである。そこから先は拡散的なので、平均的に
は 
に比例して距離が広がっていくということになる。
もっとも、これは依然として「粒子の位置がどこかわからなくなる」というこ とを意味することに変わりはない。拡散的になってしまうというのは、要する に本当にランダムになってしまうということである。もちろん、原理的には決 定論的なハミルトン系なのでランダム性なんてものはどこにもないが、しかし、 非常に近くにある2粒子の距離が、まずは指数関数的に、それからはランダム・ ウォーク的に広がるわけである。繰り返しになるが、これは物理的な系のもと もと持っている性質であり、数値計算で有限の精度でやるからそうなるとかい う問題ではない。
実際、系が熱力学的に振舞うことを保証しているのは、このように粒子の軌道 が拡散的に広がっていくという性質そのものであるはずである。近くの軌道が 指数関数的に広がる速さは熱力学的な振舞いに直接つながるわけではない。
では結局数値計算で正しいものは何で、 それはどうすれば正しいということがわかるのであろうか？ 正しいものはなにかというと、結局、個々の 粒子の位置には左右されないような統計的な性質であるというのがまあ優等生 的な答ではある。具体的には、前のほうで述べた重力熱力学的カタストロフィー や自己相似的進化がそのような系の統計的振舞いであると我々が理解できるも のの例ということになる。これらは、等質量の質点系で球対称という極めて単 純な系の場合に起きることであり、これらの条件を外した系ではもっといろい ろなことが起きるわけだが、そういうものについてもそれが統計的な性質であ る限りは「正しい」ということになる。
とはいえ、なにが統計的な性質でなにがそうでないかというのはそんなに自明 ではない。この一つの理由は、実際の系でも系を構成している粒子の数がそれ ほど多くはないことであり、もう一つは、例えば自己相似的な進化自体の結果 として粒子数が小さなサブシステムが勝手にできてきて、それがどうなるかが 系の進化を決めることになるからである。
例えば、以下のような例を考えてみよう。前に重力熱力学的振動について述べ たが、ここで、収縮および膨張自体をドライブするのは重力熱力学的不安定性 であり、これら自体はもちろん熱力学的な性質から出てくるものである。これ に対し、収縮から膨張へのスイッチングを起こすのは連星系からのエネルギー 供給である。
連星系が効率的に形成されるためには、3体衝突が起きる必要があり、そのた めには系の速度分散で正規化した密度が非常に大きい必要がある。この ためには、コアの粒子数が小さくないといけない。というか、そもそもこのような系 の粒子数に対する依存性が連星系の生成確率にあるので、中心コアが小さくな ると連星系からのエネルギー供給が無視できなくなってくるわけである。 ここで、実際に連星系ができるようになってくるのはコアの粒子数が 100 以 下になってきたところであり、同時に存在する連星は数個程度である。このた めに、連星系からのエネルギー供給には非常に大きな揺らぎが存在する。
ここで、エネルギー供給に揺らぎがあっても膨張・収縮は熱力学的に起こるの
であるから、多少の揺らぎはそのなかで吸収されると想像するかもしれない。
しかし、話はそんなに甘くはなくて、実はこの重力熱力学的振動が、エネルギー
供給が決定論的、つまり、密度の2乗といった形で与えられ、さらに粒子数
無限大の連続極限をとった場合でも、安定な
リミットサイクルを持つものではなくカオス的になるということが、フォッカー
プランク計算に人為的にエネルギー入力を入れたシミュレーションで示されて
いる
。つまり、多体計算の場合には、もともと決定論的カオスであるようなも
のにさらに粒子数が小さいことによる揺らぎが入っているわけである。このよ
うな時に、どのようにして計算結果が正しいといえるかというのは難しい。
例えば、計算精度を変えて同じような答がでるかというわけだが、粒子数が小 さいための揺らぎが系全体の振舞いを変えるので、例えば振動の位相といった ものは計算精度を変えれば全く変わってくる。多体系はカオス的で一つの粒子 の軌道には意味がないという話をしたが、このように一つの粒子の軌道が系全 体の振舞いに影響を与えるような場合にはそれでは何が正しいといえるのかと いうことが問題なわけである。
もっとも、これはもともと自由度が小さくてカオス的な系と同じ問題であり、 特に多体系の固有の問題というわけではない。結局、この場合にはさらにアン サンブル平均なり時間平均の振舞いが問題であるということになるが、カオス 的であるとそういうものがちゃんとあるかどうかも自明ではない。
まあ、そういうわけで、なにが正しくて、どうすれば正しいといえるかという
のは、結局のところ良くわからない。普通どうやっているかというと、例えば
エネルギー保存を見てシミュレーションの 
とかそんな程度で保存
していればまあ良かろうとするわけである。これが 
とかいうとまあ
計算結果を信用する気が起きない。では、 1% とか 0.1% ではいかんのかと
いうと、これはなかなかなんともいえない。例えば、重力熱力学的カタストロ
フィーによる自己相似的な中心密度の上昇を再現するだけなら、エネルギーの
誤差が例えば 0.1% と大きくてもあまり問題ではないようである。
実際問題として、エネルギーを見るというのは、特に、積分スキーム自体がエ ネルギー保存を保証しないような場合には、かなり厳しい条件になっていて大 抵うまくいく。もっとも、後述のシンプレクティック公式や時間対称型公式の ように原理的にはエネルギーの誤差が積もっていかない方法の場合には、エネ ルギー保存をみていていいのかどうかというのは明らかではない。
だいぶ話が拡散したが、ここでのもともとの議論はどうやって時間刻みを決め るかということであった。
いずれにしても軌道の時間積分の誤差をみてそれを評価したい。実用上の観点 からすると、なんらかの意味で精度を与えるパラメータが相対的な形になって いる、つまり、例えば単位系をとり直しても同じパラメータが使えることが望 ましい。これは、極めて実際的な問題としては、系の大きさや粒子数を変える たびに精度パラメータを考え直さないといけないのは面倒であり、間違いのも とになるからである。また、1つのシミュレーションの中で、系が全体として 膨張していくような場合に、それに応じて精度パラメータを変更するのは繁雑 に過ぎる。
一方、ガリレイ変換に対して不変ではないという明らかな問題から、位置や速 度の値自体を相対誤差の指標として用いることはできない。重力ポテンシャル は原理的には利用できるが、例えば銀河団の中の銀河の中にある粒子といった ものを考えると重力ポテンシャルは銀河団全体が作っているものが大きいので、 銀河団の中のどこに銀河がいるかで時間刻みが変わるという妙なことになって しまう。
そこで、重力加速度と、その高次の時間導関数だけを使って構成され、時間の 次元を持つような量が誤差の指標としては望ましい。予測子と修正子の差の絶 対値を加速度の絶対値で割ったもの
は一応そのような性質を持っている。ただし、たまたま加速度が0というよう な場合には妙なことになるので、実際に使う時には注意が必要である。
式(96)はそのような、加速度が 0 の場合でも大丈夫な形に なっている。その代わり、明示的に誤差を見ているわけではない。これは、結 局、積分スキームの誤差を評価してそれを一定に保つというよりは、むしろ単 に高階導関数の値を見て、補間多項式に対する寄与があまり大きくなり過ぎな いように適当に時間刻みを小さくするというようなものである。
例えばエネルギー保存を見た時に、この式は比較的うまく働き、例えば予測子 と修正子の差を誤差とするようなやり方とほぼ同等な性能を出す。つまり、時 間ステップの数が同じ時の誤差が同程度である。
もっとも、積分公式として4次以上のものを使った場合には、当然ではあるが 式(96)ではうまくいかない。ここでうまくいかないとは、パ ラメータを変えてタイムステップを小さくした時に、積分公式の次数から期待 されるような割合でエネルギーの誤差が小さくならないということである。 この場合にはもちろんより高次 の導関数まで入ったものを使うべきであり、そのためには式 (97)のような、誤差を明示的に評価するものをつかうのが容易 であることはいうまでもない。
「特異点」とかいうとなんかかっこいい話みたいに聞こえるが、ここではなん
ということはなくて重力ポテンシャルは距離の逆数つまり 
で、
で発散するというだけのことである。この発散するポテン
シャルがさらに引力であり、しかも古典力学系であるという事実のために、質
点系のシミュレーションでは2粒子間のポテンシャルが実際に発散し得る。こ
れに対する対策には基本的に2つある。一つは、多少物理に影響するかもしれ
ないがポテンシャルをいじって発散をなくしてしまうというものである。もう
一つは、真正直に 
ポテンシャルを使い、数値計算法を工夫することで
対応しようというものである。以下、順に解説しよう。
これは、つまり、 r が小さいところで重力の形を変えてしまって、発散が 起きないようにしようというものである。問題は、こういうことをすると系の 進化を決める物理 が変わってしまわないかということであるが、実は粒子数N がある程度大き いところではあまり物理は変わらない。
ここで、物理といっているのは2体緩和である。2体緩和は、前に説明したよう
に近接散乱よりも主に遠くの粒子との相互作用が効いている。特に、距離
程度以下、つまり、散乱角が 90度程度、あるいはそれ以上になるよう
な近接散乱は、2体緩和にはほとんど寄与していない。従って、その辺の散乱
を精密に追いかけなくても、系の進化にはほとんど影響はないと考えられる。
また、問題によってはそもそも2体緩和自体正しく表現する必要はない、ある いはむしろ積極的に押えられるものならば押えたいことがある。これは、例え ば銀河系のシミュレーションのように、現実の系では N が非常に大きくて2 体緩和のタイムスケールが非常に長いのに、シミュレーションではそれほと大 きな粒子数が使えないので2体緩和が起こってしまうというような場合である。
これらの場合には、rが小さいところで発散しないようにポテンシャルをな まらせることができる。なまらせる方法はいろいろあるが、もっとも広く使わ れているのは
というふうにする方法である。なお、文献ではこれはプラマー (Plummer)ソフ トニングと呼ばれることが多い。これは、プラマーモデルと呼ばれる自己重力 多体系の力学平衡分布モデルのポテンシャルがこの形になっているからである。 このモデルは、銀河や星団のなかの星の分布のモデルとしてそんなにいいわけ ではないが、分布関数や密度分布等が単純な初等関数で書けるという特色があ る。そのため、それほど正確な形が必要というわけではない時には銀河や星団 のモデルとして良く使われる。
歴史的には銀河団の進化のシミュレーションで、銀河をこの形のポテンシャル で相互作用する剛体粒子（といっても、重力だけで相互作用するので、相互に かさなりあうことができる）で表現した[Aar63]のがこの形が使わ れた最初であるらしい。といっても、銀河団のシミュレーションをやろうと思っ てこの形を使ったのか、この形を使った計算しかできなかったからやった計算 を「銀河団のシミュレーション」と称したのかというとどうも後者のように思 われる。
と、余談はさておき、式(98)の形のソフトニングのいいと
ころは、非常にプログラミングが楽で、計算コストも低いことである。つまり、
あらかじめ 
を適当な変数にもっていれば、加算一つだけでソフ
トニングが実現できる。さらに、非常に滑らかな関
数でソフトニングしたあとは全領域で無限回微分可能であり、高次の時間積分
法等を使っても問題が起きないということもある。
もちろん、滑らかであるということは欠点と背中合わせでもある。つまり、厳
密な 
ポテンシャルからのずれが r を大きくしてもゆっくりとしか小
さくならないのである。式(98)からすぐにわかるように、
ではポテンシャルや重力の相対誤差は 
の程
度であり、 r の2乗でしか落ちない。これはなんとなく嬉しくないような気
がするわけで、もっと高次のベキで落ちるとか、あるいは有限距離で誤差が 0
になるようなソフトニングのやりかたがいろいろ提案されている。
その中で、代表的なものはスプラインソフトニング[HK89]
と呼ばれるものである。といっても、ポテンシャル自体をスプライン関数で表
現するわけではなく、粒子の半径方向の密度分布をスプライン曲線で表現され
るとして、それを積分した形のポテンシャルを使う。ある半径 2h の外側で
は密度が 0 になるように作るので、そこでは厳密な
ポテンシャルと一致
することになる。
といっても、 r の2乗でしか落ちないことが本当に問題かどうかはあまり自
明な問題ではない。例えば実際にエネルギーの誤差を計算すれば、これは、通
常のプラマーソフトニングでは
のオーダーになっている。ではス
プラインソフトニングではどうなるかというと、これも実は 
のオーダー
になって大した違いはない。ただし、係数はスプラインのほうが多少小さい。
係数が小さいということは、誤差を同じにするなら h の値を 
の値に比べて多少大きくできるということであり、これは後で述べるようにタ
イムステップを多少大きくできる可能性があるということでもある。が、スプ
ラインの場合には関数の連続性に制限がつくわけであり、これが時間積分にど
のように影響するかはよくわかっていない。このため、どのようなソフトニン
グが「最良」であるかはまだ未解決の問題である。
なお、ソフトニングがある程度大きければ、タイムステップに実質的に下限が
できるということには注意する必要がある。今、ソフトニングが 
より
も十分に大きいとしよう。この場合、2体の近接散乱で、粒子が曲がる角度は
非常に小さくなる。また、散乱が起こっているあいだも速度はほとんど一定と
みなせる。
この時明らかなことは、タイムステップを 
の数分の1程度にとっ
てもよいということである。特にタイムステップが一定ならばこれで大きな問
題はない。ここで v は典型的な粒子間の相対速度である。
さて、これよりもずっとタイムステップを大きくとってもあまり問 題はおきないということが実験的には知られている。つまり、ある程度粒子数 が大きく、またソフトニングも大きければ、系のポテンシャルは軌道の数値 積分の精度という観点から見る限りにおいて実効的になめらかになる。このた めに時間ステップはその滑らかな軌道を積分するために必要な程度でよく、ソ フトニングの大きさとは無関係になるわけである。
もうちょっと正確にいうと、実際に滑らかになるわけではないが、タイムステッ プが大きいために高周波成分をちゃんと積分できてない（いわゆるエリアシン グエラーが入っている）分が無視できる程度に小さくなるということである。 無視できるための条件は、結局曲がる角度が小さいということである。
もちろん、この場合、エリアシングエラーによって高周波成分の寄与は増幅さ れていることには注意する必要がある。この増幅のために、各粒子のエネルギー は2体緩和で変わる分以上に変化している。これが計算結果に影響していない かどうかは確認しておく必要がある。
さらに、このようにタイムステップを大きくとれるのは、タイムステップを一 定にとった場合だけであることにも注意しておく。タイムステップを適応的に 変える場合には、結局高周波成分を拾ってそれに応じてタイムステップが変化 する必要がある。そうなっていないと、そもそもまともにタイムステップを変 化させられないからである。このために、タイムステップは通常ソフトニング サイズを分解できる程度に小さくならざるをえない。
このために、比較的大きなソフトニングを使った計算では時間刻みを下手に変 えるよりは一定にしたほうが結局良い結果がえられる。この場合に時間刻みを どのようにとればいいかもまた未解決の問題である。
ソフトニングのところでは、 
ポテンシャルのために起きる発散をそも
そも起きないようにするという話をしたが、どんな問題でもそういう扱いをし
ていいというわけではもちろんない。例えば、すでに述べたように理想化され
た質点系では重力熱力学的カタスロフィーの結果自己相似的に中心密度が増加
する。3体衝突で連星ができることによってこの自己相似的な成長は止められ
るわけだが、ソフトニングがあるとそもそも連星ができなくなってしまう。と
いうのは、連星が回りにエネルギーを有効に与えるためには、その「温度」が
回りの「温度」よりも高い必要がある、つまり、大雑把にいって連星内の軌道
運動の速度が回りの星の速度分散より大きくないといけない。これは、軌道半
径が 
よりも小さいということとほぼ同一になる。従って、ソフトニン
グがこれよりも大きいと連星ができないことになる。
このような場合には、どうしても2粒子が非常に近付いた時でもうまく計算で きるような方法が必要になる。ポテンシャル、加速度が発散するので、これは 一見原理的に不可能なような気もするが、実用的にはなんとかなる。というの は、2体が十分近付けば、その2体の相対位置については2体問題の解析解なり、 あるいはもっと高い精度が必要であれば回りの粒子からの摂動を入れた解を使 えばいいからである。近付いた粒子はいずれ離れるので近い時だけ解析解を使 えばポテンシャルが発散しても問題ではない。
もちろん、このような扱いをするためにはまずその2体を認識する必要がある が、これは単にある適当に設定した限界距離以下になるという条件でいいので それほど面倒な話ではない。
つまり、具体的には以下のようなことをする手続き群が必要になる。
さて、このセクションの表題であったところの、正則化という話はではどこに いってしまったのかというわけだが、これは以下のような話になる。
正則化の基本的な考えは、座標と時間を適当に変換して数値解の性
質をよくしようというものである。実際には3次元の運動だが、ここでは2次元
の例で説明する。2次元座標
を、
という関係を満たす新しい変
数
に変換し、さらに独立変数を時間tから rds = dt
という変数変換で得られる変数sに置き換える。すると2体問題の運動方程式
が単振動の方程式
に変換される。ここでEはシステムの全エネルギーである。単振動の方程式 の方が2体問題よりも数値的な性質がよく、少ないステップ数（大きなタイム ステップ）で精度よく積分できる。摂動も変換して扱えるので、この方法では 軌道離心率が1に近い、あるいは厳密に1であるような場合でも摂動がある解を 数値的に求めることができ、摂動が無視できるようなところは解析解でつなぐ といった余計なプログラミング上の手間は必要でなくなる。
2次元の場合にこの変換を提案したのは、 Levi Civita[Lev06] である。これは単純で2次元の場合には極めて有効であるが、例えば星団の中 の連星の軌道に使うには2次元では不便である。 3 次元の場合に拡張すること に初めて成功したのは Kustaanheimo and Stiefel[KS65] であり、 Levi Civita から半世紀 以上立った後のことである。2次元の場合の本質は、位置ベクトルを複素数と みなしてその平方根を とるということであるが、3次元ベクトルで平方根というのはうまくいかない。 Kustaanheimo and Stiefel はスピノルを使って4次元に軌道を拡張することに より、実質的に平方根演算になるような変換を定義することに成功した。
4次元で平方根というと、いわゆるハミルトンの4元数が使えないかと思う人も いるかも知れない。これはかなり最近まで良くわかっていなかったというか、 あまり誰も真剣に考えていなかっただけかも知れないが、とにかくそういう状 況であったが、最近になって可能であるということが示されている。
Kustaanheimo と Stiefel のもともとのものと4元数を使った定式化のどちら の場合でも、3次元が4次元になったぶん自由度が冗長になっているが、これは 適当に処理して問題ない。
このような正則化を行なうことで、ポテンシャルが発散していても計算精度と しては問題なく扱えるようになる。「計算精度としては」と断ったのは、実は これでは計算時間の方に問題が発生するからである。2つの粒子が遠くからた またま近付いて、双曲軌道を描いてまた離れていくようなケースでは、上の扱 いで十分である。問題なのは、連星である。
繰り返し述べたように、連星は星団の中心核の密度が上がり、3体衝突が起き るようになってくると形成される。さらに、一度形成されるとこれは他の粒子 と衝突した時に、内部の軌道運動から相手および自分の重心運動にエネルギー を与える傾向がある。
なぜ一方的にエネルギーを与えることになるかというと、連星はエネルギーを 与えた結果軌道長半径が小さくなり、その結果内部運動の速度がいっそう大き くなるからである。このこと自体が自己重力系の熱力学的な不安定性の一つの あらわれではある。
軌道長半径が小さくなるということはもちろん軌道周期が短くなるということ
である。どの程度短くなるかというと、おおよそ他の粒子が星団内を動くタイ
ムスケールの 
程度まで小さくなる。この程度で止まるのは、この
あたりで次に他の星と衝突した時に連星の重心運動がもらうエネルギーが、星
団を脱出できる程度まで大きくなるからである。
軌道周期が 
ということは、もちろん時間刻みがその程度短
くないといけないということである。従って、独立時間刻みでやったとして他の粒子全
部をあわせたのの 1000 倍の計算時間がかかるということになり、そういうも
のがごくわずかの時間でもできれば計算は全く進まなくなる。
さらに、このような非常にコンパクトな連星ができると、これに他の星が近付
いて共鳴状態になったものの計算ではタイムスケールが短いだけでなく、座標
の桁落ちが無視できなくなる。
実際には、いくつかの細かい工夫の積み重ねでこれらの問題を回避している。
まず、十分コンパクトな連星で、回りの粒子からの摂動が無視できるほど小さ い場合には、要するに単に無視すればよいわけである。この時は、連星の内部 運動はケプラー運動であり、解析解がわかっているので全く計算する必要がな くなる。
ただし、どの程度摂動が小さくなれば無視してよいかというのは理論的には微 妙な問題がある。というのは、軌道要素によって変化の摂動の大きさへの依存 性がちがうからである。もっとも鈍感なのはエネルギー（軌道長半径）であり、 これは断熱不変量なので摂動が小さくなる（変化が遅くなる）と指数関数的に 変化が小さくなる。
質点系で、連星が主に熱源として重要であるという場合には、熱源としての働 き（散乱断面積）にはエネルギー以外の軌道要素、つまり離心率やルンゲ・レ ンツベクトル（重心から近点に向いたベクトル）の影響はあまり大きくないの でエネルギーさえちゃんと計算できればあとはどうでもいい。これに対して、 現実的な効果、例えば星が有限の大きさを持ち、潮汐力によって変形したり振 動モードが励起されたり、あるいは重力波放射したりすることを考えにいれた いとすると、離心率が「正しく」計算されていることは極めて重要である。し かし、離心率の変化は摂動のベキでしか小さくならないので、これをきちんと 扱うためにはかなり摂動が小さいところまで無視できないことになる。
他の軌道要素まで正しくしようとするとさらに大変な話になるが、そこまでや
る意味はあまりないと普通は考えている。というのは、ここでの「正しく」と
いうのは所詮統計的な意味でしかないし、これは
体計算自体がそうで
あるからである。従って、ルンゲ・レンツベクトルや軌道の位相を正しく計算
しても、回りの星の軌道がそれに見合う精度では計算できていないのであまり
意味がないと考えられる。
さて、全く摂動が無視できるほど回りの星が離れてはいないが、といってまと もに計算するのは大変過ぎるという領域では、摂動論的な扱いをしたく なる。もちろん、どうやるかは天体力学の教科書をみれば書いてあるわけで、 そういうものを実装すればいいが、これは結構面倒臭い。というのは、摂動項 を求めるには連星の軌道に沿って外力を時間平均して、それから各軌道要素の 時間変化率を求めるという手続きが必要になるからである。さらに、軌道要素 とデカルト座標なり正則化座標系との相互変換をする関数も必要になる。
実効的に摂動論的計算をするのに、このような面倒な手間を一切かけないです ます画期的に安直な方法が、 Mikkola と Aarseth[MA96] によって提案された。この方法の基本的な考えは、「摂動論的扱いができる範 囲内で連星の軌道周期を遅くする」というものである。摂動項は、実時間での 軌道要素の変化率が同じになるようにスケールして運動方程式に入れる。これ によって、摂動論の範囲内では摂動論と同じ結果を、軌道を直接数値積分して 得ることができる。これはコロンブスの卵的であるが、極めて有効な方法であ る。
もちろん、こういった方法ではうまくいかない場合というものが実は存在する。 一つは、階層的な3重星、つまり連星の外側をもう一つ星が回っているような ものである。この場合には、もちろん安定であるとわかっていれば計算しなく ていいわけだが、安定性条件は数値実験によっていくつかの場合について求め られているだけで、 N 体計算で形成されるさまざまな場合について機械的 に安定性を判定する手法は存在しない。また、上に述べたようなタイムスケー ルと摂動項のスケーリングは、3重星の場合には使えない。さらに、数値的 に求まっている安定条件にしたところで、あくまでも「その時間まで計算した 範囲では壊れなかった」というだけで、もっと長時間について安定かどうかは わからない。これは、太陽系が安定かどうかというラプラスの昔からあってま だ完全には決着がついていない問題と本質的には同じ話である。
ちょっと話題がそれるが、太陽系の安定性についての最近の研究について触れ ておこう。そのへんの本を見るとラプラスが太陽系の安定性を証明したとか書 いてあるかもしれないが、現代的な理解ではこれは正しくない。この辺の詳し いことはウソを書いてはいけないので省略するが、最近の研究のハイライトは MIT のグループが摂動展開をしないで直接に外側 5 惑星の軌道の数値積分を 行ない、リヤプノフ指数が正であることを発見した[SW88] ということである。驚くべきことに、リヤプノフ指数はその逆数がわずか 2千 万年と、太陽系の年齢であるおよそ 50 億年に比べると極めて短いものであっ た。
実際の太陽系の年齢がもっと短いというのはありそうな話ではないので、リヤ プノフ指数から出てくる不安定のタイムスケールがこれほど短いというのは当 初はなかなか信用されなかった。信用されなかった理由の一つは、MIT グルー プの使っていた計算法がちょっと怪しいもので、精度が十分であるかどうかはっ きりしなかったことと、専用計算機[ADG+85]を使った大計算 で追試が困難であったことである。
90年代に入ると、汎用計算機の高速化と次節で述べるシンプレクティック法な ど数値計算法の発展の結果、追試が可能になってきた。国立天文台の中井・木 下らの追試[KN96]では、確かにリヤプノフ指数は大きいものの、太陽系は50億年とい う長期に渡って安定であるという結果になっている。つまり、リヤプ ノフ指数が単純には系の不安定性の時間尺度を与えていないのである。さらに 最近の結果では、500億年に渡って安定という結果も得られている。現実問題 としては太陽の寿命は100億年程度なので500億年という数字自体に意味がある わけではないが、いずれにしてもリヤプノフ指数から決まるタイムスケールよ りもはるかに長い時間に渡って安定であるという結果になっている。
と、そういうわけで、階層的な3重星の扱いには未解決というか解決があるか どうかもわからない問題がある。これは、数値計算法の問題というよりは、む しろカオス力学の問題であるといえよう。つまり、最終的に安定かどうかがわ からないと、どのように扱うべきかが決まらないのである。
さて、究極的に安定かどうかはともかく、階層的な3重星を安定に数値積分し ようと思えば3重星の内部運動を内側の連星と外側の連星にわける、つまりツ リー構造的な扱いをしたい。このような扱いをした場合、各自由度の運動方程 式や、摂動項がどのようになるかは、再帰的な形に書き下すことで比較的簡潔 に表現できる。[MT98,PMHM01]再帰的な扱いを可 能にしておくことで、3重星に限らず2つの連星からなる連星や、もっと粒子数 が多くて複雑なシステムも自動的に扱えるようになる。このように階層的にし た場合には、各自由度は摂動つきのケプラー運動になるのでKSによる正則化を するなりあるいは古典的な方法で積分するなり、どんな方法でも使うことがで きる。
ややこしそうな気がするのは、どのような条件でツリー構造を新しく作ったり 壊したりするかということだが、これは局所的に粒子1つだけをとったりつけ たりすることで、ある瞬間をみたときに常に最適な構造になっているとは限ら ないが、数値的に問題を起こさない程度にはよい構造を実現できることがこれ までの実績からはわかっている。
階層的でない状態にある3重星や、4個以上の系には、従来 General N-body
relaxation [Heg74] や、その簡略版である Chain
regularization[MA93b] と呼ばれる方法が使われてきた。
General N-body relaxation は、簡単にいってしまうとn個の粒子からな
る系を考えた時に、 
個のすべてのペアの相対距離に対して全部
KS 正則化をかけるような方法である。このため、方程式系は非常に冗長度が
高く、計算コストが増えるが、そのかわり座標系の切替えをまったく必要とし
ないで任意の2体の衝突を正則化できる。このために、切替えによる丸め誤差
の累積を避けることができるという潜在的な利点があることになる。
筆者個人の意見としては、しかし、どのような方法を使ったところである程度 のヒステリシスを設けることで、座標系の切替の回数は実用上問題にならない 程度には小さくできるので、計算コストが増えることを補えるだけのメリット がこれらの方法にあるとは思えない。もしも切替えによる計算精度の低下が本 当に問題になるなら、4倍精度演算を使う等して精度が落ちないようにするの は容易だからである。
さて、ここまで、数値積分自体に使う公式としては基本的にかなり古典的な PC法だけを扱ってきた。PC法は、すでに紹介したように古い歴史を持つ計算法 であり、実際その歴史は電子計算機の発明以前に遡るものである。もちろん近 年になっていろいろな発展がなかったわけではないが、我々が使っている方法 は実際に 19世紀末にアダムス等が定式化したものから大きく変わっているわ けではない。アダムス法にせよ、それを2階の方程式に特化したシュテルマー・カウ エル法にせよ、もとの微分方程式の局所的な近似解を基本的にはテイラー級数 展開のようなものに基づいて構成するというだけのものであり、それ以上に元 の系の力学的な性質を利用してよい結果を得ようというものではない。このよ うな公式をハミルトン力学系である重力多体系に適用したときには、系のエネ ルギーのような保存量には通常時間に比例する誤差が現れる。
周期解をもつ調和振動子やケプラー問題ではエネルギーのような保存量に時間 に比例する誤差が現れるのはまあ当然のことであろう。誤差が小さい、線形の 範囲であれば、最初の1周期に出た誤差と同じような誤差が繰り返しでるから である。多体系の場合には各粒子の軌道はカオス的であるとはいえ、平均的に は他の粒子全体が作る滑らかなポテンシャルのなかで準周期的な運動をするし、 時々ある他の粒子との双曲的な接近遭遇では、積分スキームや時間刻みのとり 方を決めればいつでも誤差の方向は同じになるので、やはり単位時間当たりの 誤差はおおむね一定になる傾向がある。
多体系の場合にはどうせカオス的だし、エネルギーが保存しないからといって どれほど悪いかというのはそれ自体がよくわからない問題であるのはすでに述 べた通りだが、前に述べた弱い摂動を受けた連星系を長時間積分するような場 合だとエネルギーや角運動量が保存しないのは嬉しくない。というのは、それ らさえ保存すればあとはどうでもいいというところがなくはないからである。 従って、そのような「良い」公式はないかというのが実用上の問題になる。
経験的には非常に古くから、低次であるにもかかわらず特別によい振舞いをす る公式があるということは、重力多体系、分子動力学計算、プラズマ物理など ハミルトン系の数値計算が使われるどの分野でも知られていた。それぞれ別の 名前がついており、また分野によっては（丸め誤差の範囲で）同じものがいく つかの名前で呼ばれているが、実体はみな同じものでいわゆるリープフロッグ 公式である。 (Verlet, 2nd order ABM, 2nd order Stömer-Cowell あたり が古くからある名前の代表的なものである)一つの形は以下のようなものであ る。
ここで、添字はステップである。 
は中間点での値というこ
とになる。
これでは速度と位置がずれた時間でしか定義されない
が、出発用公式として
を使い、さらに終了用公式として
を使うことで最初と最後を合わせることが出来る。この形は、実は
と数学的に等価である。こう書くと PC法みたいに見える。また、以下のようにも書ける
これは一見なんだか良くわからない。
さらにまた、速度を消去して
, 
, 
の関係式の形で
書いてあることもあるかもしれない。
この公式は局所誤差が 
、大域誤差が
である。
つまり、普通にいう2次精度である。
局所誤差という観点からはこれは決して良い公式というわけではないが、
現実にはこの公式は非常に広く使われている。
これは、別にもっとよい方法を知らないからとかではなく、実はこの方法がい くつかの意味で非常に良い方法であるからである。いくつかの意味とは、例え ばエネルギーや角運動量のような保存量が非常によく保存するということであ る。これらの保存量については多くの場合に誤差がある程度以上増えない。
この、ある程度以上誤差が増えないというのはめざましい性質である。通常の方 法では、エネルギーの誤差は時間に比例して増えていく。従って、長時間計算 をしようとすればそれだけ正確な計算をする必要がある。ところが、エネルギー の誤差は溜っていかないのならば、かならずしも精度を上げる必要はないとも 考えられる。
もちろん、エネルギーが保存していればそれだけで計算が正しいということに はならない。が、理論的には、いくつかの重要な結果が得られている。まず、
シンプレクティック公式については、以下のようなことが知られている
それぞれについての詳しい解説は、それだけでかなりの量になるので「数理科 学」の吉田による解説[Yos95]等をみていただくことにして、ここ ではそれぞれの考え方に基づいてより高次の公式を構成するいくつかの方法に ついて簡単にまとめる。
とはいえ、能書きを聞いていても良くわからないので、実例を見てみよう。まず、簡単な 例として調和振動子
を リープフロッグ と 2階のルンゲクッタで解いた例を示す。初期条件は 
で、時間刻みは 
である。
軌道とエネルギーを図に示す。非常に特徴的なのは、 リープフロッグ ではエネル ギーが周期的にしか変化しないのに対し、ルンゲクッタでは単調に増えている ことである。
図 7: 調和振動子の数値積分。左：軌道。右：エネルギー。破線は2次のル
ンゲクッタ、実線はリープフロッグの結果
ルンゲクッタでは単調に変化するというのは前に説明した通りである。エネル ギーが増えるか減るかは自明ではなく、ルンゲクッタでもどんな公式を使って いるかや時間刻みの大きさによるが、2次の場合には公式や時間刻みによらず に必ずエネルギーが増える方向にいく。 これに対し、リープフロッグではエネルギーが変化していない。これはど ういうことなのであろう。
実は、この調和振動子の場合には、リープフロッグ公式は以下の量
を保存するということを確かめることが出来る。つまり、 
で与えら
れる位相平面の上で考えると、リープフロッグ公式の解は上の式で与えられる楕円
の上にのっているのである。このために、エネルギーの誤差がある値よりも大
きくなり得ないことになる。
では、非線形振動ではどうだろう？簡単な例として
をリープフロッグと 2階のルンゲクッタで解いた例を示す。初期条件は 
で、時間刻みは 
である。
軌道とエネルギーを図に示す。調和振動の場合と同様に、リープフロッグではエネル ギーが周期的にしか変化しないのに対し、ルンゲクッタでは単調に増えている。
図 8: 非線形振動の数値積分。左：軌道。右：エネルギー。破線は2次のル
ンゲクッタ、実線はリープフロッグの結果
リープフロッグ 公式は、すでに述べたようにハミルトン系に対してエネルギー誤 差が有界に留まるという大きな特長がある。が、2次精度でしかない。もっと 高次の方法はないのだろうか、またあるとすればどういう原理で作れるのだろ うか。
高次の方法を構成する一つのアプローチが、シンプレクティック公式といわれ るものである。これはなにかというと、積分公式がシンプレクティック写像に なるように作るということである。
シンプレクティック写像とは、要するに正準変換のことである。なぜそれで はカノニカル公式といわないでシンプレクティック公式という難しげな言葉を 使うのかは知らないので聞かないでほしい。正準変換とはなにかというのは、 まあ直観的には力学系を不変に保つ、つまり変換前の座標系で求めた軌道を 変換したものと、変換後の座標系での力学系の軌道が同じになるようなもので ある。
ハミルトン力学系の解そのもの（ある時刻 t での座標から、 t+h での座 標に移す変換）もシンプレクティックである。まあ、だから、シンプレクティッ クになっているような積分公式は、そうでないものに比べてなんとなく力学系 の性質にあっているような気はする。
なお、リープフロッグがシンプレクティックであることを示すのは、これが ある瞬間の加速度で速度だけを進める写像と、ある瞬間の速度で位置だけを進 める写像であることに注意すると、それらがそれぞれシンプレクティックであ ることを示せばいいことになって結構簡単である。さらに、このことから、リー プフロッグのようにこの2つを組み合わせて構成される写像、特にリープフロッ グを刻み幅を変えて繰り返し適用する写像はすべてシンプレクティックであ ることになる。
で、いいたかったことは何かっていうと、上の リープフロッグ公式はこのシ ンプレクティック性を満たしているということだった。以下、高次の公式にはどんなものがあるかという話をしておく。
陽解法の組み立て方は、基本的にはルンゲ・クッタ・ナイストレム型の公式、 すなわち、2階の方程式に特化したルンゲ・クッタ型の公式で、係数をシン プレクティック性を満たすように決めるということである。これはここ10年で 無数に論文がでた。4次、あるいは6次の公式としては、吉田や鈴木による作用 素分解に基づく公式が良く知られていて、広く使われていると思う。
これらの方法の原理は、要するに先に述べたようにリープフロッグをタイム ステップを変えていくつか組合わせるというものである。うまくタイムステッ プを組み合わせると誤差の高次の項を消すことができるわけである。7段6次や 15段8次の公式等が吉田によって導かれている。鈴木による再帰的な公式は、 原理的にいくらでの高次のものが構成できるという意味では美しいが、指数関 数的に段数が増えるので実用性からはちょっと疑問な気もする。
なお、実は陽解法はハミルトニアンが 
の形の場合にしか使え
ないが、大抵の問題はこう書けることはいうまでもないであろう。
また、 RKN 系の公式を力任せに構成する試みもあり、4次から8次までの公式 が作られている。計算精度という観点からはこちらのほうがリープフロッグを組合 わせるものよりも何桁も良いようである[SSC94]。
さて、シンプレクティックであるということと、「良い」ということの間には
ちゃんとした理論的な関係があるのであろうか？一応あるということになって
いる。つまり、あるハミルトニアン H で表される系をシンプレクティック
な p 次の公式を使って積分した解は、別のハミルトニアン 
で与えられ
るシステムの厳密解になっていて、H と 
の間に
という関係がある（
を求める数列が収束すれば）ということがわかってい
る。
なお、例えば上の調和振動子の場合には実際に数列が収束し、
が求まっ
ている。が、これは極めて例外的な場合で、一般の場合には収束するかどうか
も明らかではないようである。
収束するかどうか明らかではないのでは、使っていいことがあるという保証は ないではないかと思うかもしれない。理論的にはその通りなのであるが、実際 にはいろいろな問題に適用されて、従来の方法よりも高い精度が得られるとい うことが確認されている。
さて、式(115) をみるとわかるように、シンプレクティック公式
に付随するハミルトニアン
には時間刻み h が入っている。従って、
h をふらふら変えると 
も変わって、結果的に求まった数値解は一つの
力学系の軌道ではない変なものになってしまう。ということは、可変時
間刻みにするとシンプレクティック公式はうまく働かないのではないかという
ことが想像される。
実はその通りで、例えばシンプレクティックで埋め込み型のルンゲクッタとい うものを作って、実際に時間刻みを変えてみた人がいる。その結果、普通のル ンゲクッタよりも良くならないということが発見された [CSS93]。多くの 問題でこれは致命的な欠陥となるので、さまざまな対応法が精力的に研究され ているのが現状である。が、あまり一般的にうまくいく方法というのは見つかっ ていないようである。つまり、万能な公式というのはシンプレクティック公式 に関する限り知られていない。
なお、シンプレクティック公式において普通に時間刻みを変えると、それはシ ンプレクティックではなくなるということが Skeel と Gear により一般に証 明されている[SG92]。彼らの証明は、大雑把にいうと時間刻みも含めて考えると、写 像が時間方向に歪んだもの（もとの位置によって時間刻みが変わるため）にな り、そのために元の公式がシンプレクティックであれば時間刻みを変えると必 ずシンプレクティックでなくなるというものである。
それでは、時間刻みを変えるのは絶対に不可能かというと、以下のような考え かたでの研究が進められている。ハミルトニアン H を
と複数の項の和にわけ、それぞれについて違う時間刻みを与えることで実効的 に時間刻みを変えるというものである。（例えば [SB94]）
上のようなハミルトニアンの書き換えは、自由度が小さい系ならできなくはな いかもしれないがここでテーマにしている重力多体系では実際的ではない。と いうわけで、現在のところ、古典的なリープフロッグ公式がソフトポテンシャル と組み合わせて一定時間刻みで使われている他は、シンプレクティック公式が そのまま重力多体系に使われることはほとんどない。
もう一つのシンプレクティック公式の問題点は、「写像」であるということか
ら一段法、具体的にはルンゲクッタ型の公式であるので、局所誤差に対する計
算量という観点からは線形多段階法に比べて必ず悪いということである。実は
線形多段階法を持ち出すまでもなく局所誤差に対して最適化したルンゲ・クッ
タ法に比べても、驚くほど誤差の係数は大きい。例えば 8 次の公式の場合、
15段の吉田による公式の誤差の係数は、2階の問題に最適化された非シンプレ
クティックな Dormand らによる9段公式（実質 8 段 8 次）の
倍にも
なる。ステップ当たりの計算量が約2倍であることを考慮すると、同じ計算時
間での局所誤差は100万倍も大きいということになる。まあ、8次の公式なので、
精度を同じにするとして計算時間は6倍程度で済むが、よほど長時間の積分で
なければ非シンプレクティックな公式のほうが実は良かったりする。4次の場
合も、同様に良い非シンプレクティックな公式と普通に使われているリープフ
ロッグ3回の組合わせでは2桁ほど局所誤差に差がある。
なお、シンプレクティックであってなおかつ局所誤差が小さい公式というもの も存在はしている。代表的なものは陰的ガウス公式である。これは、ガウス・ ルジャンドル型の数値積分公式から導かれるいわゆる collocation 法 といわ れるクラスの陰的ルンゲクッタ公式である。ルンゲ・クッタなので積分区間の 途中の何点かで関数の値を計算する。陰的ガウス公式ではこの点をルジャンド ル多項式の 0 点にとり、計算した導関数値から補間多項式をつくってそれを 区間で積分することで次のステップでの解を求める。さらに、collocation 法 であるとは、途中の点で導関数値を計算するのにその点での近似解がいるが、 その解に補間多項式を積分したものを使うことを指している。このためにこの 方法では補間多項式の係数についての陰的な方程式を解く必要がある。
この方法の重要な利点は、段数に対して2倍の次数を達成できることである。 つまり、4点を使う公式で8次が実現できる。陽的な方法ではシンプレクティッ クでなくても8次には9段が最低限必要であり、それよりもだいぶ少ない。ま た、局所誤差の係数も極めて小さい。
しかしながら、陰的公式であり何らかの方法で方程式を解かないといけない。 自由度が大きい時には何らかの方法といっても直接代入以外は現実的ではない ので、まあまあよい初期推定をもとにすることで数回で反復を止められるように しないと結局かえって遅くなることにはなる。とはいえ、原理的にはいくらで も高次の方法が機械的に構成でき、さらに計算コストも少なくなる可能性があ るという意味で魅力的な方法である。単なる2体ケプラー問題で 8 次のシンプ レクティック公式同士を比べた時には、直接代入で収束させたガウス公式がもっ ともよいという結果が得られている[SSC94]。
この方法のさらにもう一つの利点は、ハミルトニアンが分離可能でなくても使 えるということである。
前節で、 シンプレクティック公式がそのまま重力多体系に使われることはほ とんどないと書いたが、これは「そのまま」でなければ使われている例 [WH91]もあるということでもある。ここではそれがどんな 例か説明しておく。とりあえず、どういうものに使いたいかというと、例によっ て摂動を受けたケプラー問題である。この方法自体は、重力多体系といっても 恒星系ではなく、太陽系の惑星の長時間積分のために開発されてきたものであ るが、すでに述べた連星の扱い等にも応用されるようになってきているし、ま た他の分野への応用もひろいと考えられるので少し詳しく見ていくことにする。
リープフロッグと、その組合せからなるようなシンプレクティック公式の基本的な 考えは、ハミルトニアン
に対して、pに対する以下の変換
と q に対する同様な変換がそれぞれシンプレクティックなので、それらを 順番に適用したものもシンプレクティックになるということであった。
さて、摂動をうけたケプラー問題では、ハミルトニアンは
と書ける。 惑星系なら
が太陽からの重力で、 
は自分以外の惑星
からの寄与ということになる。このようにハミルトニアンを分解した時に、
の解はそれ自体シンプレクティックであり、また 
だけを考えた
というマッピングもシンプレクティックである。従っ
て、この2つを組み合わせて積分公式をつくることができる。ここでやってい
ることは、結局普通のリープフロッグでは 
と直
線で動かす代わりに、ポテンシャル 
に沿って動かすというだけである。
ケプラー問題は解析的に解けるので、このようなやり方で太陽中心の軌道を極
めて高精度に積分できることになる。 
の寄与については普通にやった
のでは2次精度になるが、普通のリープフロッグと全く同様にステップサイズを変
えたものを組み合わせることで高次にすることもできる。
実用的にこの方法を使うためには、ケプラー問題を高速に解く方法が必要にな る。これは、実空間でまともに解く方法と、 KS 変換した座標系で解く方法の それぞれで極めて高速に解ける方法が提案されている。太陽系の惑星のように 軌道離心率がそれほど大きくない場合には実空間で扱うのが簡単でいい。離心 率が大きい場合にはKS 変換の利用も検討はしてみるべきであろう。
なお、このような「ハミルトニアンの分割」という見栄えのいい方法ではない が、太陽重力とそれ以外の惑星からの摂動を別に扱おうという考え方は天体力 学では非常に古くからあり、エンケの方法として知られている。古くからある だけになにをもってエンケの方法というかも本によって微妙に違うので、以下 には最近の榎森ら[EIN93]による定式化の簡略版を述べる。
もともとの運動方程式を
という形に書く。ここで
は太陽からの重力であり、
は外力である。外力のない2体問題の解析
解が
であたえられる時、
は以下の
方程式に従うことになる
この方程式の右辺はルンゲ・クッタやエルミート法ではそのまま数値積分でき
る。 線形多段階法の場合は少し複雑になるが、うまく使う方法はある。この
方法では、数値積分する
の大きさは、もともとの
より
もずっと小さいので、それだけ計算精度があがる。特に、外力が非
常に小さい時にはタイムステップを長くとれる。なお、2次の MVS はこの方程式を
リープフロッグで積分していることに相当する。
こちらの定式化の利点は、MVSに比べて高次で局所誤差の小さい方法が使える ことである。そのかわり、摂動項に対してはシンプレクティックでなくなって いるので、時間に比例する誤差が完全にはなくならない。
シンプレクティック公式は、とりあえずうまく動くのでいろいろな分野で使わ れるようになりつつあるが、欠点もいくつかあるのはすでにみた通りである。 大雑把にいえば、シンプレクティック公式の利点はほぼそのままにした上で、 欠点に対応できそうなのが対称型公式ということになる。
まず対称型公式とはなにかということだが、まずここでいう時間反転に対する 対称性というものを定義しておこう。時間反転に対する対称性とは、常微分方 程式
に対する一段法
が、
を満たすこと、つまり、直観的には、ある微分方程式系があって、それを1ス テップ分数値的に軌道を進めたとする。で、そこから逆に戻ると厳密に元の値に戻る ということである。（ここでは丸め誤差はないものとする）
で、一段法の場合には、この意味で対称なものを対称型公式ということにする。
具体例で考えてみよう。例えば前進オイラー法は対称型ではない。これは、いっ た先で導関数を計算すれば、一般にもとのところでの導関数とは違うから当然 であろう。前進オイラー法が対称ではないのだから、その逆写像である後退オ イラー法も対称型ではないことになる。
では、対称型にはどういうものがあるかということだが、明らかに対称型であ るものとしては、台形公式
がある。これが上の対称性を満たしていることはいうまでもない。
さて、なぜこの型の公式を考えるかということであるが、少なくとも実験的に は以下のようなことが知られている。[HMM95]
前のほうの性質はシンプレクティック公式と同じであるが、後の方の性質はシ ンプレクティック公式よりもある意味でよいものである。
さて、一階の方程式に対する一段法という制限をつけた場合、つまり、ルンゲ クッタ法の場合、対称な公式はかならず陰的公式になる。逆に、陰的公式であ れば対称型のものを作るのは容易である。具体的には、前に出てきた陰的ガウ ス法は対称型公式でもある。
一階の系に対する通常の一段法では、これよりよい方法は多分ない。それでは、
以下、順番に考えていくことにする。
実は、すでに述べたように、リープフロッグ公式は対称型でありしかも陽公 式である。で、リープフロッグの組合わせで作られる公式も、実は すべて対称型になっている。
さて、一段法という枠を外してみよう。
まず、線形多段階法に対する対称性の定義を与える必 要がある。ここでは、以下のように定義しておこう。 時間刻みが一定の時、2階の方程式用の線形多段階法は、一般に
という形に書くことができる。 
なら陽的であるし、そうでな
ければ陰的公式である。ここで、対称型であるとは、係数が
を満たすということである。リープフロッグをこの形にしたものは、
というものであり、これが上の対称性の定義を満たしていることはいうまでも ないであろう。
上の意味で対称な公式は、実際に時間反転に対して対称であるということに注
意しておこう。つまり、 
から上の式で 
を求めたとする。時間を逆転させた解を考えて、
から 
を求めるということを考えても、使
う式が時間が逆転していないものと同じになっているのである。
このような公式が存在することは 1970 年代から知られていたが、注目される ようになったのはシンプレクティック法の研究が盛んになった 1990 年代に入っ てからである。特にトロント大の Quinlan と Tremaine[QT90]が14次までの対称型 公式を導いている。
なお、非常に最近になって、国立天文台の福島によって、Quinlan & Tremaine によるものよりもはるかに良い性質を持つものが導かれている。[Fuk98]
ただし、これらの公式に共通する欠陥として、線形解析では安定であっても、 非線形問題（例えば、単純なケプラー問題）に適用すると共鳴により不安定な 振動を起こすことがあるという困難があり、これはまだ未解決の問題である。[Fuk99,ET99] このために6次以上の公式は現在のところ実用になっていない。
福島はこれを解決する方法として、「超陰的方法」というものを提案している。
理屈はある意味簡単で、
を求めるのにもっと未来の情報まで使うと
いうものである。もちろん、このためには複数ステップを「一度に」解く必要
があり、初期値問題をあたかも境界値問題であるかのように連立方程式として
解くことになる。こういった方法が実用になるかどうかはまだこれからの研究
の発展を見ないとわからない。
さて、前に紹介したエルミート公式だが、一般の次数の場合にはそうではない
が、広く使われている4次の公式は実は時間対称である。ちなみにシンプレク
ティックではない。といっても、PEC モードで使っていたのでは陰公式が厳
密に満たされてないので、時間対称公式のもつ良い性質ははっきりとは現れな
い。極限修正モード、あるいは 
や
モード
で使うことで実質上時間対称性を実現することができる。
とはいえ、普通に独立時間刻みと合わせて使っていたのでは、シンプレクティッ
クの場合と同様に時間刻みを変えることが時間対称性を破る。また他の粒子の
軌道を予測する時には過去の情報しか使っていないので、ここでも対称性が破
れている。このために、通常のエルミート公式と独立時間刻みを合わせた
やりかたでは、
にしても特にいいことはない。
いいことがあるのは、可変時間ステップとはいってもほとんどの場合に時間刻
みは一定であるような場合である。具体的には、惑星形成のシミュレーション
では、多数の微惑星が重力相互作用しながら太陽のまわりをほとんど円軌道で
回っている。これらの時間刻みはほとんど一定であり、そのために
にするだけで劇的に計算精度が向上する。さらに、誤差のほとんど
が太陽重力から来ているので、粒子間の相互作用は一度計算するだけで、太陽
重力についてだけ反復することで精度をあげることができる。[KYM98]
さて、シンプレクティック公式は時間刻みの変更と相性が悪いのはすでに述べた とおりだが、対称型公式も普通にやったのでは同じである。が、対称型公式で は「時間刻み対称化」によって良い性質を取り戻すことができることがわかっ ている。ここでは時間刻み対称化について簡単にまとめる。
可変刻みの場合、対称型公式であったとしても普通は「時間を反転して逆に積 分すれば元に戻る」という性質を満たさなくなる。理由は簡単で、時間刻みが 同じになるとは限らないからである。つまり、時間刻みは普通は前のステップ での誤差とかそういったものを使っているので、行きと帰りでは違ったものに なる。
時間刻みが対称でないというのが永年誤差がでる原因だとすれば、対称にして
しまえば永年誤差がなくなるかも知れない。対称にするとは、時刻
で計算した時間刻み 
と 逆に戻る時に
で計算した
時間刻み 
が等しくなるようにするということである。このため
には、時間刻みが時刻 
と時刻 
の情報を同じ重みで使っていれば
良い。一つの可能性は以下のようなものである。
ここで 
は時刻 
での位置と速度であり、 
は適当な時
間刻みを与える関数である。つまり、積分の始点と終点で時間刻みをそれぞれ
計算して、その平均を取るわけである。こうすれば、逆に戻った時にも時間刻
みが等しくなることが保証される。ここで注意しないといけないのは、
は 
の関数であるということである。このために、式
(131) は 
を陰的に決める方程
式になっていて、何らかの方法で解く必要がある。これには疑似ニュートン法
などいろいろな方法が使える。
Cano と Sanz-Serna [CSS97]は、周期運動するハミルトン系について、(a)時 間刻み一定のシンプレクテック公式, (b)時間ステップ対称化した対称型公式、 (c)ハミルトニアンを保存するように構成された公式の3つについて、任意の量 の誤差が最大でも時 間の1次でしか増大しないということを証明した。これは時間刻み対称化にあ る程度の理論的裏付けを与えるものと解釈できよう。
自己重力系を応用先として考える時には、時間対称化が有効なのはまずは連星 系等である。共鳴状態の3体のような、直接に数値積分する以外に特に有効な 手段がなく、しかも時間刻みを変えることは絶対的に必要な場合によい結果を 与えるからである。しかし、独立時間刻みと組み合わせることはできないだろ うか？
Quinn et al. [QKSL97]は、リープフロッグは対称化した上で独立時間刻みを使う
ことが可能であることを示している。これは、基本的にはリープフロッグが陽公式
であり、位置は予測子が修正されないでそのまま使われるためである。なお、
速度は修正子によって修正される。従って、対称性を保証するためには、時間
刻みが速度ではなく位置だけに依存するようになっている必要がある。彼らは
時間刻みを 
に比例させることで速度に依存しない時間刻みを
作っている。が、彼らの実験結果を見る限りでは、やはり密度だけによって時
間刻みを決めることには無理があるのか、非常によい結果が得られているとは
いい難い。そのせいかどうか知らないが論文もまだ通っていないようである。
ここまでは時間積分の手法の話をしてきたが、ここでは加速度の計算の話にな る。ツリーなどの近似算法を使ってある時刻での加速度を高速に計算するとい うのは次章の話題として、ここでは独立時間刻みをもう少し進めて、ある粒子 への力をいくつかの成分にわけ、それぞれに違う時間刻みを与えるという方法 について説明する。 「いくつか」といっても、実用になっているのうは近くの粒子からの力と遠く の粒子からの力の2つに分けるもので、理論的にはそれ以上分けてもあまりメ リットがないということがわかっている。[MH88]
分ける原理は、粒子毎に適当な半径 r を考えて、ある時刻でそのなかにい
た粒子はネイバー粒子であるとし、リストをつくって憶えておくというもので
ある。その時刻では、全粒子からの力をネイバーからの力 
と残りからの
力 
に分けて計算しておく。タイムステップはこの2つのそれぞれについ
て計算する。普通は
のほうが短いステップサイズになるので、時間が
のタイムステップを超えるまでは 
はステップごとに計算するが、
の方は外挿する。超える手前の時刻で、
を計算しなおして修正子
を適用すると同時に、ネイバー粒子のリストも作り直す。
なお、実際にこのアルゴリズムを線形多段階法と組み合わせるのはかなり面倒
である。というのは、
を新しく計算した時点でネイバー粒子のリストも
変わるので、単純な差分によって予測子を構成できないからである。広く使わ
れていたAarseth によるプログラム[Aar85]では、入れ替わった分
の補正を入れ替わった粒子について解析的に高次の導関数を計算して行なう。
この処理はかなり繁雑なものになる。
但し、4次のエルミート法の場合には予測子に前のステップの情報を使わない ので、このような補正は不要になりプログラムを単純化できている。[MA92]
このやりかたでどの程度得になるかというと、理論的および数値的な解析の結
果では、ネイバー数は
程度が最良で、計算量の比は
程度になる。[MH88,MA92]
粒子数の 
乗というとたいしていいことがないように思うかもしれないが、
これは必ずしもそうではない。というのは、コストにつく係数が小さいために
比較的粒子数が小さいところ (100程度)でも既に効果があるからである。粒子
数が
程度では、計算コストに 10倍程度の差が出てくる。
通常の計算機を使った高精度の重力多体計算では、現在現実的に可能な程度の 粒子数ではこのネイバースキームを使った方法がもっとも高速であるようであ る。
ここまでの話では、どういう計算機の上で実行するかということ はあまり考えないで、精度を確保し、計算量（浮動小数点演算の数）を減らす にはどういうことを考えればいいかということを見てきた。 30 年前であれば ここまでで話はほとんどおしまいであり、あとはかなり細かいプログラミング テクニックの問題ということになった。
しかし、最近30年間の計算機の発展は、ある意味ものごとをより面倒にする方 向に来ている。というのは、20年前であれば高速な計算機というのはベクトル プロセッサであり、アルゴリズムがベクトル化できるかどうかで性能が数十倍 違ったし、現在ではもっとも高速な汎用計算機は分散メモリ型並列計算機になっ ていてここでは並列化できるかどうかで性能が千倍以上違ってくるからである。
独立時間刻みを使うような問題の並列化、特に分散メモリの計算機での並列化 については、しかし、現在のところあまり研究がないのが現状である。これが どうしてかは良くわからないが、次に述べるツリー法に比べて本質的に難しい というのが一つの理由であろう。まず、計算量が多い割に扱う粒子数が小さく、 効率的な並列化のためにはプロセッサ間の通信や同期などが非常に高速である 必要がある。現在、並列計算機の主流になりつつある分散メモリ型並列計算機 や PC クラスタがもっとも苦手とするタイプの計算なのである。これに対し、 古典的なベクトル計算機や共有メモリ型の並列計算機ではあまり難しいことを しなくても性能がでる。
ベクトル計算機や共有メモリ型の並列計算機で 単純な独立時間刻みを使う場合、力の計算には全粒子からの力を受けるという高い並列 性があるので、そこのところを並列化すれば十分である。これに対してネイバー スキームとなってくると並列度が下がるので、もう少し増やしたい。 これにはブロックステップ法というものが使われる。[McM86]こ の方法では、基本的に時間刻みを2のベキ乗に量子化することで、なるべく沢 山の粒子が同じ時刻をとるようにする。同じ時刻にくるものは並列に計算でき る。
分散メモリの場合、今のところどのように並列化すべきかということがそもそ も良くわかっていない。提案されているのは例えば全プロセッサが全粒子のコ ピーを持つ方法である。[SB99]今のところ、この方法 ではあまり高い性能が出ていない。使っているプラットホームは Cray T3D/E であり、プロセッサ間通信が非常に速いものであるので、それでも性能が出て いない理由はよくわからない。理論的には、この方法で通信が遅い PC クラス タで良い性能を出すのはきわめて難しい。
汎用並列計算機での効率的な実装の研究が進んでいないもう一つの理由は、最 後に述べる専用計算機が単純な独立時間刻みに対しては極めてよい成果を出し ているということかもしれない。これについては後でもう少し詳しく議論する。
