太陽系の数値モデルとは、一連の数式であり、これを解くことで、惑星のおおよその位置を時間の関数として求めることができる。このようなモデルを作成しようとする試みが、より広範な天体力学という分野を確立した。このシミュレーションの結果は、過去の測定値と比較して精度を確認し、将来の位置を予測するために利用できる。したがって、その主な用途は天体暦の作成である。
シミュレーションは、直交座標系または球面座標系のどちらでも実行できます。前者はより簡単ですが、計算負荷が非常に高く、電子計算機でしか実用的ではありません。そのため、以前は後者のみが使用されていました。厳密に言えば、後者も計算負荷がそれほど低くはありませんでしたが、簡単な近似から始めて、必要な精度に達するまで摂動を加えることが可能でした。
本質的に、この太陽系の数学的シミュレーションは、N体問題の一種です。記号Nは天体の数を表し、太陽、8つの惑星、数十の衛星、無数の小惑星、彗星などを含めると、かなり大きくなります。しかし、太陽が他の天体に及ぼす影響は非常に大きく、他のすべての天体が互いに及ぼす影響は非常に小さいため、この問題は解析的に解ける2体問題に還元できます。各惑星の結果は軌道であり、その位置を時間の関数として単純に記述したものです。これが解かれると、衛星と惑星が互いに及ぼす影響が小さな補正として加えられます。これらは完全な惑星軌道に比べると小さいものです。一部の補正は数度になる場合もありますが、測定は1秒角よりも高い精度で行うことができます。
この方法はシミュレーションにはもはや用いられていませんが、比較的単純な主解を用い、いくつかの大きな摂動を加えることで、比較的容易に目的の惑星位置にたどり着けるため、近似的な天体暦を求めるには依然として有用です。ただし、摂動理論は非常に高度な数学を必要とするという欠点があります。
現代的な手法は、3次元空間における数値積分から成ります。まず、関係する各物体の位置 ( x、y、z ) と速度 ( v x、v y、v z ) の高精度な値から始めます。各物体の質量も分かっている場合は、ニュートンの万有引力の法則から加速度 ( a x、a y、a z ) を計算できます。各物体は互いに引き合い、全体の加速度はこれらの引力の合計になります。次に、小さな時間ステップ Δ tを選択し、ニュートンの運動の第 2 法則を適用します。加速度に Δ tを掛けると速度の補正値が得られ、速度に Δ tを掛けると位置の補正値が得られます。この手順を他のすべての物体に対して繰り返します。
その結果、すべての物体の位置と速度の新しい値が得られます。次に、これらの新しい値を使用して、次の時間ステップ Δt について計算全体を最初からやり直します。この手順を十分に繰り返すことで、すべての物体の位置の時間変化を記述したデータが得られます。
この方法の利点は、コンピュータにとって非常に簡単な作業であり、摂動を決定するための複雑で困難な手順を省き、すべての天体について同時に非常に正確な結果が得られることです。欠点は、最初に非常に正確な数値を用意しなければ、結果が時間とともに現実から乖離してしまうこと、x、y、z座標が得られるため、使用する前に、より実用的な黄道座標または赤道座標に変換する必要がある場合が多いこと、そして、すべてを計算するか、何も計算しないかのどちらかしかないことです。ある特定の時刻におけるある惑星の位置を知りたい場合、他のすべての惑星とすべての中間的な時間ステップも計算する必要があります。
前節では、加速度が小さな時間ステップΔtの間一定であると仮定し、計算は単にV×ΔtをRに加えるだけで済むと考えました。しかし実際には、Δtを非常に小さくしてステップ数が膨大になる場合を除き、これは当てはまりません。なぜなら、位置は常に加速度によって変化しますが、加速度の値は瞬間的な位置によって決まるからです。明らかに、完全な積分が必要となります。
いくつかの方法があります。まず、必要な方程式に注目してください。
この方程式は、特定の物体jに対して、1 から N まで移動するすべての物体i の加速度を表しています。これはベクトル方程式なので、X、Y、Z の各成分について 3 つの方程式に分割し、次の式を得ます。
追加の関係性
、
YとZについても同様です。
前者の式(重力)は難しそうに見えるかもしれませんが、計算は問題ありません。後者の式(運動法則)はより単純に見えますが、計算はできません。コンピュータは積分ができず、微小な値を扱うことができないため、dt の代わりに Δt を使用し、結果として得られる変数を左辺に移動させます。
、 そして:
aは依然として時間の関数であることに注意してください。これらを解く最も簡単な方法は、オイラーのアルゴリズムです。これは本質的には、上述の線形加算です。一般的なコンピュータ言語で1次元のみに限定すると次のようになります。
a.old = gravitationfunction(x.old) x.new = x.old + v.old * dt v.new = v.old + a.old * dt
本質的に、タイムステップの全期間で使用される加速度は、タイムステップの開始時の加速度と同じであるため、この単純な方法では高い精度が得られません。開始値と期待される(摂動のない)終了値の平均である平均加速度を取ることで、はるかに優れた結果が得られます。
a.old = gravitationfunction(x.old) x.expect = x.old + v.old * dt a.expect = gravitationfunction(x.expect) v.new = v.old + (a.old + a.expect) * 0.5 * dt x.new = x.old + (v.new + v.old) * 0.5 * dt
もちろん、中間値を用いることでさらに良い結果が得られる場合もあります。これは、ルンゲ・クッタ法を用いる場合によく見られる現象で、特に4次または5次のルンゲ・クッタ法が最も有効です。最も一般的に用いられる方法は、長期的なエネルギー保存性に優れていることから、リープフロッグ法です。
全く異なる方法として、テイラー級数を用いる方法があります。その場合、次のように記述します。
しかし、r に関してのみ高次の導関数まで展開するのではなく、次のように記述することで r と v (つまり r') に関して展開することができる。そして、因数fとgを数列で書き出します。
加速度を計算するには、各物体が他の物体に及ぼす重力による引力を考慮に入れる必要があります。その結果、シミュレーションにおける計算量は物体の数の二乗に比例して増加します。物体の数を2倍にすると、計算量は4倍になります。シミュレーションの精度を高めるには、小数点以下の桁数を増やすだけでなく、時間ステップを小さくする必要がありますが、これもまた計算量を急速に増加させます。当然ながら、計算量を減らすための工夫が必要となります。ここでは、そうした工夫のいくつかを紹介します。
最も重要なコツは、既に上で述べたように、適切な積分方法を用いることです。
単位の選択は重要です。SI単位系では、極端に小さい値と極端に大きい値が出てしまうため、すべての単位を1に近い値にスケーリングする必要があります。例えば、太陽系内の距離を表すには、天文単位が最も分かりやすいでしょう。これを怠ると、計算中に浮動小数点オーバーフローやアンダーフローが発生し、シミュレーションが途中で中断される可能性が非常に高くなります。たとえオーバーフローが発生しなくても、切り捨て誤差によって精度が損なわれる可能性があります。
Nが大きい場合(太陽系シミュレーションではそれほどではないが、銀河シミュレーションではより顕著)、動的な天体群を作成するのが一般的である。その時点で計算されている基準天体から特定の方向かつ遠距離にあるすべての天体をまとめて、それらの重力による引力をグループ全体で平均化する。
閉鎖系の全エネルギーと角運動量は保存量です。これらの量を各タイムステップ後に計算することで、変化が小さい場合はステップサイズΔtを大きくし、変化し始めた場合は小さくするようにシミュレーションをプログラムできます。また、前述のように物体をグループ分けし、近い物体よりも遠い物体に大きく、つまり少ないタイムステップを適用することも可能です。
特定の物体が基準物体に近づくと加速度が過度に急激に変化するため、小さなパラメータeを導入するのが一般的です。
最高精度が求められる場合、計算ははるかに複雑になります。彗星の場合、放射圧やガス抵抗などの非重力的な力も考慮する必要があります。水星やその他の惑星の場合、長期計算では相対論的効果を無視することはできません。また、全エネルギーももはや一定ではありません(線運動量を持つ4元ベクトルエネルギーは一定であるため)。光速が有限であるため、古典的および相対論的な光時間効果も考慮する必要があります。惑星はもはや粒子として扱うことはできず、その形状と密度も考慮する必要があります。例えば、地球の扁平化は歳差運動を引き起こし、それが自転軸の傾きの変化を引き起こし、すべての惑星の長期的な動きに影響を与えます。太陽系の安定性の欠如のため、数千万年を超える長期モデルは不可能です。