N体シミュレーション

バーンズ・ハット・シミュレーションなどのツリー法では、通常、オクツリーを使用してボリュームを立方体のセルに分割し、近くのセルの粒子間の相互作用のみを個別に処理する必要があります。遠く離れたセルの粒子は、遠く離れたセルの重心を中心とする単一の大きな粒子（または低次の多重極展開）としてまとめて処理できます。これにより、計算する必要のある粒子ペアの相互作用の数を大幅に削減できます。粒子間の相互作用の計算によってシミュレーションが過負荷にならないようにするには、セルあたり多数の粒子を含むシミュレーションの密集部分で、セルをより小さなセルに細分化する必要があります。粒子が均等に分布していないシミュレーションでは、 CallahanとKosarajuの十分に分離されたペア分解法により、固定次元で反復ごとに最適なO( n log n )時間が得られます。  

もう一つの可能​​性は、粒子メッシュ法である。この方法では、空間をメッシュ上に離散化し、重力ポテンシャルを計算する目的で、粒子はメッシュの周囲の2×2の頂点に分割されていると仮定する。ポテンシャルエネルギーΦは、ポアソン方程式を用いて求めることができる。

$${\displaystyle \nabla ^{2}\Phi =4\pi G{\rho },\,}$$

ここでGはニュートン定数であり、${\displaystyle \rho }$は密度（メッシュ点における粒子の数）です。高速フーリエ変換は、ポアソン方程式が単純な形式となる周波数領域に移行することで、これを効率的に解くことができます。

$${\displaystyle {\hat {\Phi }}=-4\pi G{\frac {\hat {\rho }}{k^{2}}},}$$

どこ${\displaystyle {\vec {k}}}$は共動波数であり、ハットはフーリエ変換を表す。${\displaystyle {\vec {g}}=-{\vec {\nabla }}\Phi }$重力場は、を掛けることで求めることができる。${\displaystyle -i{\vec {k}}}$そして逆フーリエ変換を計算する（または逆変換を計算してから他の方法を用いる）。この方法はメッシュサイズによって制限されるため、実際にはより小さなメッシュ、あるいは他の手法（ツリーや単純な粒子間アルゴリズムとの組み合わせなど）を用いて小規模な力を計算する。場合によっては適応メッシュが用いられることもあり、その場合、シミュレーションの密度の高い領域ではメッシュセルがはるかに小さくなる。

太陽系内の天体の軌道をかなり正確に推定するために、いくつかの異なる重力摂動アルゴリズムが使用されている。

人工衛星を凍結軌道に乗せることはよくあります。地球を周回する人工衛星の軌道は、地球の中心を周回する2体楕円軌道を基本として、地球の扁平率、太陽と月の引力、大気抵抗などによる微調整を加えることで正確にモデル化できます。そのため、人工衛星の実際の軌道を計算することなく、凍結軌道を見つけることが可能です。

小型惑星、彗星、あるいは長距離宇宙船の軌道は、太陽の周りを回る2体楕円軌道を起点とし、既知の軌道にあるより大きな惑星の重力による微調整を加えることで、多くの場合正確にモデル化できる。

粒子系の長期的な経路の特性の一部は直接計算できる。個々の粒子の実際の経路を中間段階として計算する必要はない。こうした特性には、リアプノフ安定性、リアプノフ時間、エルゴード理論に基づく様々な測定値などが含まれる。

ソフトニングは、N体シミュレーションにおいて、粒子が別の粒子に近づきすぎた場合（力が無限大になる場合）に数値発散を防ぐために用いられる数値的なトリックです。これは、各粒子の正則化された重力ポテンシャルを次のように修正することによって得られます。

$${\displaystyle \Phi =-{\frac {1}{\sqrt {r^{2}+\epsilon ^{2}}}},}$$

（1/rではなく）${\displaystyle \epsilon }$は軟化パラメータです。シミュレーションが現実的なものとなるよう、軟化パラメータの値は十分に小さく設定する必要があります。

N体シミュレーションの最も単純な実装では、${\textstyle n\geq 3}$これは、軌道上の物体の単純な伝播です。単純とは、軌道上の物体に作用する力は、物体同士が互いに及ぼし合う重力のみであることを意味します。C ++などのオブジェクト指向プログラミング言語では、伝播に必要な基本的な数学的構造とデータコンテナを確立するために、いくつかの定型コードが役立ちます。具体的には、状態ベクトル、したがってベクトル、そしてこのデータと軌道上の物体の質量を含む基本的なオブジェクトです。この方法は、他のタイプの N 体シミュレーションにも適用できます。電荷を持つ点質量のシミュレーションでは、同様の方法を使用しますが、力は電場の相互作用による引力または斥力によるものです。いずれにせよ、粒子の加速は、力ベクトルの合計を粒子の質量で割った結果です。

$${\displaystyle {\vec {a}}={\frac {1}{m}}\sum {\vec {F}}}$$

粒子の運動学的データを格納するための、プログラム的に安定かつスケーラブルな方法の一例として、固定長配列の使用が挙げられます。最適化されたコードでは、メモリの割り当てが容易になり、消費されるリソースを予測できます。以下のC++コードをご覧ください。

struct Vector3{double e [ 3 ] = { 0 };Vector3 () {}~ Vector3 () {}inline Vector3 ( double e0 , double e1 , double e2 ){this -> e [ 0 ] = e0 ;this -> e [ 1 ] = e1 ;this -> e [ 2 ] = e2 ;}};struct OrbitalEntity{double e [ 7 ] = { 0 };OrbitalEntity () {}~ OrbitalEntity () {}inline OrbitalEntity ( double e0 , double e1 , double e2 , double e3 , double e4 , double e5 , double e6 ){this -> e [ 0 ] = e0 ;this -> e [ 1 ] = e1 ;this -> e [ 2 ] = e2 ;this -> e [ 3 ] = e3 ;this -> e [ 4 ] = e4 ;this -> e [ 5 ] = e5 ;this -> e [ 6 ] = e6 ;}};

OrbitalEntity状態ベクトルを格納するのに十分なスペースがあることに注意してください。ここで、

• ${\textstyle e_{0}=x}$、デカルト空間におけるオブジェクトの位置ベクトルの投影${\displaystyle \left[1\;0\;0\right]}$
• ${\textstyle e_{1}=y}$、デカルト空間におけるオブジェクトの位置ベクトルの投影${\displaystyle \left[0\;1\;0\right]}$
• ${\textstyle e_{2}=z}$、デカルト空間におけるオブジェクトの位置ベクトルの投影${\displaystyle \left[0\;0\;1\right]}$
• ${\textstyle e_{3}={\dot {x}}}$物体の速度ベクトルをデカルト空間に投影したもの${\displaystyle \left[1\;0\;0\right]}$
• ${\textstyle e_{4}={\dot {y}}}$物体の速度ベクトルをデカルト空間に投影したもの${\displaystyle \left[0\;1\;0\right]}$
• ${\textstyle e_{5}={\dot {z}}}$物体の速度ベクトルをデカルト空間に投影したもの${\displaystyle \left[0\;0\;1\right]}$

さらに、OrbitalEntity質量値のための十分なスペースがあります。

一般的に、N体シミュレーションは、何らかの運動方程式に基づくシステムになります。これらのシステムのほとんどは、シミュレーションを「シード」するための初期構成に依存します。重力や電位に依存するシステムなどでは、シミュレーションエンティティにかかる力はその速度に依存しません。したがって、シミュレーションの力のシードには、初期位置のみが必要ですが、これだけでは伝播はできません。初期速度が必要です。恒星の周りを公転する惑星を考えてみましょう。惑星は運動していませんが、主星の重力の影響を受けています。時間が経過し、時間ステップが追加されると、惑星は加速度に応じて速度を増します。ある瞬間には、$t_n$ただし、隣接する質量による物体の合力加速度は、その速度に依存しないが、時間ステップでは$t_n+1$その結果生じる位置の変化は、伝播が速度に本質的に依存するため、大きく異なります。以下で使用するシンプレクティックオイラー法などの基本的な伝播メカニズムでは、$t_n+1$速度のみに依存する$t_n$位置のずれは次のように計算されます

$${\displaystyle {\vec {r}}_{t_{n+1}}={\vec {r}}_{t_{n}}+{\vec {v}}_{t_{n}}\cdot \Delta t}$$

加速なしで、${\textstyle {\vec {v}}_{t_{n}}}$静止しているが、位置だけを見ている観測者の視点からすると、速度の変化を見るには2つの時間ステップかかる。

太陽系のようなシミュレーションは、中心星から惑星に相当する点質量までの平均距離を求めることで実現できます。コードを簡潔にするため、長半径と平均速度に基づく非厳密なアプローチを採用します。これらの天体を構成する前に、メモリ領域を確保する必要があります。拡張性を確保するため、mallocコマンドを使用できます。

OrbitalEntity * orbital_entities = malloc ( sizeof ( OrbitalEntity ) * ( 9 + N_ASTEROIDS ));orbital_entities [ 0 ] = { 0.0 , 0.0 , 0.0 , 0.0 , 0.0 , 0.0 , 1.989e30 }; // 太陽に似た恒星orbital_entities [ 1 ] = { 57.909e9 , 0.0 , 0.0 , 0.0 , 47.36e3 , 0.0 , 0.33011e24 }; // 水星に似た惑星orbital_entities [ 2 ] = { 108.209e9 , 0.0 , 0.0 , 0.0 , 35.02e3 , 0.0 , 4.8675e24 }; // 金星に似た惑星orbital_entities [ 3 ] = { 149.596e9 , 0.0 , 0.0 , 0.0 , 29.78e3 , 0.0 , 5.9724e24 }; // 地球に似た惑星orbital_entities [ 4 ] = { 227.923e9 , 0.0 , 0.0 , 0.0 , 24.07e3 , 0.0 , 0.64171e24 }; // 火星に似た惑星orbital_entities [ 5 ] = { 778.570e9 , 0.0 , 0.0 , 0.0 , 13e3 , 0.0 , 1898.19e24 }; // 木星に似た惑星orbital_entities [ 6 ] = { 1433.529e9 , 0.0 , 0.0 , 0.0 , 9.68e3 , 0.0 , 568.34e24 }; // 土星に似た惑星orbital_entities [ 7 ] = { 2872.463e9 , 0.0 , 0.0 , 0.0 , 6.80e3 , 0.0 , 86.813e24 }; // 天王星に似た惑星orbital_entities [ 8 ] = { 4495.060e9 , 0.0 , 0.0 , 0.0 , 5.43e3 , 0.0 , 102.413e24 }; // 海王星に似た惑星

ここでN_ASTEROIDS、 は一時的に 0 のままになる変数ですが、ユーザーの裁量で将来的に多数の小惑星を含めることができます。シミュレーションの設定における重要なステップは、シミュレーションの時間範囲を設定することです。${\displaystyle t_{0}}$に${\displaystyle t_{\text{end}}}$増分時間ステップも同様${\displaystyle dt}$これによりシミュレーションが進行します。

double t_0 = 0 ;double t = t_0 ;double dt = 86400 ;double t_end = 86400 * 365 * 10 ; // 約10年分（秒単位）double BIG_G = 6.67e-11 ; // 重力定数

上記で確立された位置と速度は、以下の条件を満たすものとして解釈される。${\displaystyle t=t_{0}}$。

シミュレーションの範囲は論理的に次の期間となる。${\displaystyle t_{0}\leq t<t_{\text{end}}}$。

シミュレーション全体は無数の時間ステップから構成される可能性があります。基本的なレベルでは、各時間ステップでは、各物体について以下を計算します。

• 身体にかかる力。
• 加速度（${\displaystyle {\vec {a}}}$);
• 速度（${\displaystyle {\vec {v}}_{n}={\vec {v}}_{n-1}+{\vec {a}}_{n}\cdot \Delta t}$);
• および位置（${\displaystyle {\vec {r}}_{n+1}={\vec {r}}_{n}+{\vec {v}}_{n}\cdot \Delta t}$)

上記は、whileループで非常に簡単に実装できます。${\displaystyle t}$前述の範囲内に存在する：

while ( t < t_end ) { for ( size_t m1_idx = 0 ; m1_idx < 9 + N_ASTEROIDS ; m1_idx ++ ) { Vector3 a_g = { 0 , 0 , 0 };for ( size_t m2_idx = 0 ; m2_idx < 9 + N_ASTEROIDS ; m2_idx ++ ) { if ( m2_idx != m1_idx ) { Vector3 r_vector ;r_vector.e [ 0 ] = orbital_entities [ m1_idx ] .e [ 0 ] - orbital_entities [ m2_idx ] .e [ 0 ] ; r_vector.e [ 1 ] = orbital_entities [ m1_idx ] .e [ 1 ] - orbital_entities [ m2_idx ] .e [ 1 ] ; r_vector.e [ 2 ] = orbital_entities [ m1_idx ] .e [ 2 ] - orbital_entities [ m2_idx ] .e [ 2 ] ;​​​double r_mag = sqrt ( r_vector . e [ 0 ] * r_vector . e [ 0 ] + r_vector . e [ 1 ] * r_vector . e [ 1 ] + r_vector . e [ 2 ] * r_vector . e [ 2 ]);倍速加速度= -1.0 * BIG_G * ( orbital_entities [ m2_idx ]. e [ 6 ]) / pow ( r_mag , 2.0 );Vector3 r_unit_vector = { r_vector . e [ 0 ] / r_mag , r_vector . e [ 1 ] / r_mag , r_vector . e [ 2 ] / r_mag };a_g.e [ 0 ] + = acceleration * r_unit_vector.e [ 0 ] ; a_g.e [ 1 ] + = acceleration * r_unit_vector.e [ 1 ] ; a_g.e [ 2 ] + = acceleration * r_unit_vector.e [ 2 ] ; } }​​​​​​orbital_entities [ m1_idx ] .e [ 3 ] + = a_g.e [ 0 ] * dt ; orbital_entities [ m1_idx ] .e [ 4 ] + = a_g.e [ 1 ] * dt ; orbital_entities [ m1_idx ] .e [ 5 ] + = a_g.e [ 2 ] * dt ; }​​​for ( size_t entity_idx = 0 ; entity_idx < 9 + N_ASTEROIDS ; entity_idx ++ ) { orbital_entities [ entity_idx ]. e [ 0 ] += orbital_entities [ entity_idx ]. e [ 3 ] * dt ; orbital_entities [ entity_idx ]. e [ 1 ] += orbital_entities [ entity_idx ]. e [ 4 ] * dt ; orbital_entities [ entity_idx ]. e [ 2 ] += orbital_entities [ entity_idx ]. e [ 5 ] * dt ; } t += dt ; }

シミュレーションにおいて内側の4つの岩石惑星に焦点を当てると、上記の伝播によって生じる軌道は以下のようになる。

• フォン・ヘルナー、セバスチャン（1960年）。 「シュテルンハウフェンのケルパー問題の数理統合。I」。Zeitschrift für Astrophysik (ドイツ語)。50：184。ビブコード：1960ZA....50....184V。
• フォン・ヘルナー、セバスチャン（1963年）。 「Sternhaufen におけるKörper問題の数理統合。II」。 Zeitschrift für Astrophysik (ドイツ語)。57:47。ビブコード: 1963ZA....57...47V。
• Aarseth, Sverre J. (2003).重力N体シミュレーション：ツールとアルゴリズム. Cambridge University Press . ISBN 978-0-521-12153-8。
• Bertschinger, Edmund (1998). 「宇宙における構造形成のシミュレーション」。Annual Review of Astronomy and Astrophysics . 36 (1): 599–654 . Bibcode : 1998ARA & A..36..599B . doi : 10.1146/annurev.astro.36.1.599 .
• ビニー、ジェームズ、トレメイン、スコット（1987）。銀河力学。プリンストン大学出版局。ISBN 978-0-691-08445-9。
• Callahan, Paul B.; Kosaraju, Sambasiva Rao (1992). "多次元点集合の分解とk近傍法およびn体ポテンシャル場への応用（予備版）". STOC '92: Proc. ACM Symp. Theory of Computing . ACM.。