N体重力サンドボックスベータ

N体重力サンドボックスについて

このサンドボックスは、画面上のすべての天体の組について万有引力を省略せずに足し合わせる(直接和。計算量は天体数の2乗に比例)。位置と速度の更新は速度ベルレ法の kick-drift-kick 形式で行う。速度を半ステップ分進め、位置を1ステップ動かし、新しい位置で力を計算し直してから残りの半ステップ分の速度を加える。

右側の「エネルギー誤差」は、開始時の全エネルギー(運動エネルギー+位置エネルギー)と現在値のずれを百分率で示す。重力だけがはたらく系では全エネルギーは一定のはずなので、この数字はそのまま数値積分の誤差になる。合体や天体の追加があると基準を取り直す。

単位はゲーム用に決めたもの(太陽級の恒星の質量が4000、重力定数の既定値が100)。軌道の形や保存則は現実と同じだが、画面上の秒や質量を天文学の単位に換算することはできない。

オイラー法だと惑星が逃げていく理由

最も単純な前進オイラー法は、今の速度で位置を、今の加速度で速度を同時に更新する。楕円軌道では1周ごとにエネルギーを少しずつ足してしまう。誤差が打ち消し合わず、同じ向きに積み重なるからだ。

計算してみると差ははっきりする。GM = 1、距離1から接線速度0.8で出発する軌道(離心率0.36、周期約3.96)を時間刻み0.01で積分した。オイラー法では1周後にエネルギーが22.9%増え、10周後に57%、100周後に89%増えた。軌道は周回ごとに膨らみ、脱出に向かう。同じ条件で kick-drift-kick は途中のどの時点でも誤差0.021%を超えず、100周しても増えなかった。

この差は刻みを細かくして埋まる種類のものではない。ベルレ系の方法はシンプレクティック積分法で、本当のエネルギーにごく近い「修正エネルギー」を厳密に保存する。そのため誤差は1周の中で上下するだけで蓄積しない(Hairer・Lubich・Wanner 2003)。Verlet が1967年に分子動力学でこの方法を使ったのも、何十万ステップにわたってエネルギーが漏れては困るからだった。このツールの「積分法」を「オイラー (比較用)」に変えると、同じ軌道で「エネルギー誤差」が増え続け、軌道が外へほどけていく様子を直接見られる。

エネルギー誤差を大きくするもの

ベルレ法でも刻みが粗ければ誤差は増える。このツールは1フレーム(1/60秒)を部分ステップに分け、天体が少ないほど細かく(最大16)分ける。「時間の速さ」を上げると1部分ステップが受け持つ時間が長くなり、誤差が大きくなる。太陽系プリセットで速さ1と4を比べてみるとよい。

最大の誤差は2天体の近接通過で生まれる。力は距離の2乗に反比例するので、すれ違う瞬間に加速度が1ステップの間に激変する。「ソフトニング」は距離を√(r² + ε²)に置き換えてこの発散を抑える。代わりに、εより近い距離ではもうニュートン重力ではない。ソフトニングを0にして恒星のすぐ横に惑星を投げると、ありえない速さで弾き飛ばされることがある。これは物理ではなく数値誤差だ。

天体が座標5万より遠くへ飛ぶと削除される。合体・破砕・追加・削除のときと同じく、天体の数が変わった時点で基準エネルギーを取り直すので、誤差表示は残った天体のエネルギーを基準に続く。

三体問題と8の字解

2天体の運動はケプラー軌道として厳密に解けるが、3つ目が加わると一般的な解の公式はない。ほとんどの初期条件で運動はカオス的で、最初の位置がわずかに違うだけで、しばらく後の軌跡がまったく違ってくる。同じクラスの恒星を3つ似た速さで投げると、しばらく絡み合ったあと1つが弾き出され、残る2つが連星になることが多い。

それでも美しい周期解は存在する。等しい質量の3天体が1本の8の字曲線の上を追いかけ合う解だ。Moore(1993)が数値計算で見つけ、Chenciner と Montgomery(2000)が変分法で存在を証明した。G = m = 1 の単位でよく引用される初期条件(Simó の計算値)は、2天体が(±0.97000436, ∓0.24308753)、中央の天体の速度が(−0.93240737, −0.86473146)、外側の2天体はそれぞれその半分の逆向き、周期は約6.3259。

この初期条件を刻み0.001で積分すると、また違いが出る。kick-drift-kick では10周期後も中央の天体が原点に戻り、エネルギー誤差は10⁻¹¹%程度だった。オイラー法では1周期で3%、10周期で20%ずれ、8の字が崩れた。「8の字軌道 (三体)」プリセットは、この初期条件をツールの単位に合わせて拡大したものだ。3つの星の質量はM = 4000、長さの単位Lは表示領域の短辺の半分の0.42倍(画面が小さいときは、最も近づいても星どうしが触れない大きさまで広げる)、速度は√(G·M/L)倍で、周期は6.3259 × L ÷ √(G·M/L)になる。積分法をオイラーにしてプリセットを読み込み直すと、8の字が崩れる過程も比べられる。

衝突モデルは近似

「現実的」な衝突では、衝突エネルギーと結合エネルギーの比で合体・かすめて通過・部分破壊・完全破壊に分かれる。質量・運動量・重心は保存するが、破片の数や大きさ、超新星・キロノヴァが起きる質量のしきい値はゲーム用の値だ。惑星形成や巨大衝突の研究で使うSPHのような流体計算とは別物である。

衝突判定は、部分ステップの間に2つの天体がまっすぐ動いたとみなし、その区間で最も近づいた距離が半径の和より小さかったかを見る(線分の最近接点検査)。そのため、とても速い小天体でも1ステップで相手を飛び越えてすり抜けることはない。すれ違った後で見つかった衝突は、質量中心はそのままに最も近づいた瞬間の配置へ戻して衝突の向きを決める。

例

使い方

1. プリセット(太陽系・連星系・原始惑星系円盤など)を選ぶか、空の画面から始める。
2. 天体の種類を選び、ドラッグして離すと投げられる。軽くタップすると最も重い天体を回る円軌道の速さで置かれる。
3. 重力定数・時間の速さ・ソフトニング・積分法(速度ベルレ/オイラー)を変え、軌道と「エネルギー誤差」の変化を見る。
4. 必要なら3Dに切り替え、軌道面が傾いた系を回して見る。

主な用語

シンプレクティック積分法

相空間の体積を保つ積分法。エネルギー誤差が蓄積せずに振動する。速度ベルレ法やリープフロッグ法が代表例。

ソフトニング(ε)

重力の計算で r の代わりに√(r² + ε²)を使うこと。近接通過で力が無限大に近づくのを防ぐ。

コレオグラフィー解

複数の天体が同じ閉曲線を時間差でたどる周期解。8の字解が最初の例として有名。

よくある質問

補足

2Dモードではすべての天体が z = 0 の平面にあり、z方向の力は0になる。3Dから2Dに切り替えると天体は平面に押しつぶされるので、軌道が変わる。