万有引力の数値シミュレーションにRK4（4次ルンゲ・クッタ）で十分すぎる理由：2階常微分方程式とテイラー展開の観点から
重力多体系や天体軌道のシミュレーションを自作したことがある人なら、一度は「ルンゲ・クッタ法（RK4）は非シンプレクティックだから、長期間回すとエネルギーが保存されずにドリフトする」という議論を目にしたことがあるはずです。
確かに理論的にはその通りです。
しかし、実際のコードに落とし込み、適切な初期値と十分な時間刻みを与えて検証してみると、
万有引力による滑らかな天体軌道に対して、RK4は実用上十分すぎるほどの精度を発揮できる
ことが分かります。
この記事では、なぜ万有引力（逆二乗法則）のシミュレーションにおいてRK4がこれほど強力に機能するのかを、数学的な構造、初期値、時間刻み、そして実際の長期積分結果から整理してみたいと思います。
1. 運動方程式の次数とRK4の近似次数
古典力学における万有引力による質点の運動は、ニュートンの運動方程式によって
$$\ddot{\mathbf{r}}=-\frac{GM}{|\mathbf{r}|^3}\mathbf{r}$$
と記述できます。
これは位置 $\mathbf{r}$ に関する2階の常微分方程式です。
ただし、RK4は2階微分方程式へ直接「4次法」として適用するわけではありません。
通常は、
$$\mathbf{v}=\dot{\mathbf r}$$
を導入して、
$$\dot{\mathbf r}=\mathbf v$$
$$\dot{\mathbf v}=-\frac{GM}{|\mathbf r|^3}\mathbf r$$
という1階の連立常微分方程式へ変換してからRK4を適用します。
状態ベクトルを
$$\mathbf y=\begin{pmatrix}\mathbf r\ ,\mathbf v\end{pmatrix}$$
とすると、
$$\frac{d\mathbf y}{dt}=\begin{pmatrix}\mathbf v\ ,-\frac{GM}{|\mathbf r|^3}\mathbf r\end{pmatrix}$$
となります。
RK4は、この1階連立常微分方程式に対して、テイラー展開の4次までの情報を利用することで、局所誤差 $O(h^5)$、大域誤差 $O(h^4)$ を持つ数値積分法です。
つまり重要なのは、
「2階の方程式だから4次のRK4が優れている」
という単純な次数の比較ではありません。
重要なのは、
万有引力による運動方程式を1階連立ODEへ変換したとき、その滑らかな時間発展をRK4が4次精度で追跡できる
という点です。
逆二乗型の重力場は、衝突や特異点の近傍などを除けば非常に滑らかな系です。
この性質とRK4の高い局所精度が組み合わさることで、天体軌道の数値計算において非常に高い精度を得ることができます。
2. 「RK4ではエネルギーが維持できない」の意味
RK4はシンプレクティック積分器ではありません。
したがって、理論上は非常に長時間積分すれば、エネルギー誤差が完全にゼロのまま維持されるわけではありません。
ここは重要です。
しかし、
「非シンプレクティックである」ことと「実用上すぐに軌道が崩れる」ことは同じではありません。
実際の数値計算では、
- 初期位置
- 初期速度
- 重力パラメータ $GM$
- 時間刻み $\Delta t$
- 計算する時間スケール
がすべて結果に影響します。
特に初期値が理論軌道と整合していない場合、RK4の精度以前に、その時点ですでに異なる軌道を計算することになります。
逆に、力学的に整合した初期値を与え、時間刻みを十分小さく取れば、RK4の数値誤差を非常に小さく抑えることができます。
したがって、
RK4の長期積分を評価するときは、「シンプレクティックかどうか」だけを見るのではなく、実際の時間刻みと初期値でどの程度の誤差が発生するのかを検証することが重要です。
3. 万有引力において適切な初期値とは？
適切にケプラー方程式を解くことができるのであれば、初期値はケプラー軌道から取得することをお勧めします。
特に重要なのが初期速度です。
軌道要素から速度の大きさを求める公式は理論的には正しいのですが、実装上では、
- 速度の大きさ
- 速度の方向
- 座標系
- 真近点角
- 離心率
- 近日点・遠日点付近での数値精度
などをすべて正しく整合させる必要があります。
そこで、ケプラー方程式から求めた軌道上の2点を利用して、その位置差から初期速度を求める方法があります。
例えば時刻 $t$ と、周期Pに対して十分小さな時間差 $\Delta t$ を持つ
$$\mathbf r(t)$$
$$\mathbf r(t+\Delta t)$$
を求め、
$$\mathbf v(t)\approx\frac{\mathbf r(t+\Delta t)-\mathbf r(t)}{\Delta t}$$
（※ $\Delta t$ を小さくしすぎると丸め誤差が支配的になるため、周期 $P$ の 10万分の1〜100万分の1程度が実装上の最適点となる。）
とすることで、軌道上の位置に整合した速度ベクトルを得ることができます。
このようにして得られた初期値は、単に速度公式から大きさだけを与える場合と比較して、実装上の座標・方向の不整合を避けやすくなります。
4. 実装の美学：無駄を削ぎ落としたRK4
天体シミュレーションにおいて、シンプレクティック法は確かに強力です。
しかし、実装の簡潔さと汎用性においてRK4は非常に優れています。
必要なのは、状態ベクトルを正しく管理することだけです。
/**
* RK4による1ステップの更新
* @param {number} t 現在時刻
* @param {Vector} y 現在の状態（位置と速度）
* @param {number} dt タイムステップ
* @param {Function} f 加速度を返す関数 (t, y) => dy/dt
*/
function rk4Step(t, y, dt, f) {
const k1 = f(t, y); // この瞬間の傾き
const k2 = f(t + dt * 0.5, y + dt * 0.5 * k1); // k1方向への半ステップ未来の傾き
const k3 = f(t + dt * 0.5, y + dt * 0.5 * k2); // k2方向への半ステップ未来の傾き
const k4 = f(t + dt, y + dt * k3); // k3方向への1ステップ未来の傾き
return y + (dt / 6.0) *
(k1 + 2.0 * k2 + 2.0 * k3 + k4);
}
ここで厳密には、k1～k4 は「加速度そのもの」ではなく、状態ベクトル $\mathbf y$ の時間微分です。
状態ベクトルを位置と速度で構成した場合、その中に
- 位置の時間微分＝速度
- 速度の時間微分＝加速度
の両方が含まれています。
この構造によって、位置と速度を一体として4次精度で更新できます。
🌍 補足：円軌道で理解するRK4のk1〜k4の意味
作図説明
K1：接線方向（現在の速度方向）＋中心方向（現在の加速度方向）
K2：K1 による予測位置での接線＋中心方向
K3：K2 による予測位置での接線＋中心方向
K4：K3 による予測位置での接線＋中心方向
RK4:(K1+2K2+2K3+K4)/6のベクトル
この図を見れば分かる通り
K1(1点の傾き)のオイラー法に比べて
RK4(4点の傾き)はより精密に1ステップ先を
予測していることが分かりますね
RK4の4つの傾き（k1〜k4）は、未来の軌道方向を4点でサンプリングして混ぜ合わせることで、曲率のある軌道を高精度に追跡するための構造です。
円軌道はもっともシンプルな「方向ベクトルが常に変化する」系なので、k1〜k4の意味を直感的に理解できます。
k1：現在状態での傾き
現在位置における速度と加速度から、
「今この瞬間の軌道がどちらへ進むか」
を評価します。
k2：k1の方向に半ステップ進んだ状態
半ステップだけ未来へ進んだ仮想状態を作り、その位置における速度と加速度を評価します。
軌道上では位置が変化するため、中心方向も変化します。
k3：k2の方向に半ステップ進んだ状態
k2と同じ時刻
$$t+\frac{\Delta t}{2}$$
を評価しますが、k1から予測した状態ではなく、k2を使って予測した状態を評価します。
そのため、k2とk3は同じ時刻であっても一般に異なる状態となります。
k4：1ステップ未来の状態
k3を利用して1ステップ先の仮想状態を作り、その地点での傾きを評価します。
RK4：k1～k4の重み付き平均
$$\frac{1}{6}(k_1+2k_2+2k_3+k_4)$$
という重み付き平均を取ることで、4つの評価点から1ステップ全体の変化を高精度に近似します。
🌌 なぜk2とk3の時刻が同じなのか？
k2とk3はどちらも
$$t+\frac{\Delta t}{2}$$
の時刻を評価します。
しかし、
- k2はk1を使って半ステップ先を予測
- k3はk2を使って半ステップ先を予測
するため、評価する状態が異なります。
つまり、
同じ時刻でも、そこへ到達すると予測した経路が異なる
ということです。
これによって、単純な「現在と未来の2点」だけでは捉えにくい軌道の曲率を、RK4は内部的に複数回評価できます。
✨ まとめ：円軌道はRK4の構造が見えやすい
円軌道では、
- 加速度が常に中心向き
- 速度方向が連続的に変化する
- 軌道が滑らか
- 状態ベクトルを直感的に理解できる
ため、RK4の内部計算を理解する教材として非常に適しています。
RK4の4つの評価点が、軌道の時間発展をどのように追跡しているのかを視覚的に確認できます。
🔬 実測：高離心率軌道を3750年間積分してみる
ここまでの説明を、実際の数値計算で検証してみます。
今回使用した軌道は、
$$a=17.81466\ {\rm AU}$$
$$e=0.95$$
$$P=75\ {\rm 年}$$
という高離心率軌道です。
特に離心率 $e=0.95$ は、かなり極端な楕円軌道です。
この軌道をRK4で50周、つまり、
$$50\times75=3750\ {\rm 年}$$
相当積分しました。
その結果、古典軌道要素を再評価すると次のようになりました。
| 軌道要素 | 理論値 | 50周後 | 再現率 |
|---|---|---|---|
| 近日点距離 | 0.8907329 AU | 0.8903233 AU | 99.954 % |
| 近日点周期 | 75.00000 年 | 74.98280 年 | 99.977 % |
| 軌道長半径 | 17.81466 AU | 17.81460 AU | 100.000 % |
| 離心率 | 0.9500000 | 0.9500228 | 100.002 % |
特に軌道長半径は、
$$17.81466\ {\rm AU}\rightarrow17.81460\ {\rm AU}$$
であり、相対誤差は約
$$-0.00034%$$
です。
離心率も、
$$0.9500000\rightarrow0.9500228$$
で、変化量は約
$$2.28\times10^{-5}$$
に収まっています。
さらに10周、20周、50周の結果を比較すると、軌道要素が周回数に比例して大きくドリフトしている様子は確認できませんでした。
これは、
RK4だから長期積分では必ず軌道が崩壊する
という単純な理解ではなく、
適切な初期値と時間刻みを選択すれば、RK4でも実用的な天体シミュレーションにおいて非常に高い長期精度を得られる
ことを示す実例です。
もちろん、この結果だけから「任意の時間スケールでRK4が必ず安定する」と結論することはできません。
しかし少なくとも、今回のような高離心率軌道について3750年相当の積分を行っても、主要な軌道要素がほぼ理論値の周辺に維持されていることは、RK4の実用性を示す十分に興味深い結果です。
🌠 では、シンプレクティック法は不要なのか？
もちろん、そういう意味ではありません。
シンプレクティック積分器には、長期間のハミルトン系において非常に優れた性質があります。
特に、
- 数百万年
- 数千万年
- 数億年
といった極端な長期積分や、エネルギー・位相空間構造の保存を重視する計算では、シンプレクティック法が有力な選択肢になります。
一方で、
「シンプレクティックでなければ天体シミュレーションには使えない」
というわけではありません。
実際のシミュレーションでは、対象とする時間スケール、必要な精度、時間刻み、計算量、実装の複雑さを総合的に判断する必要があります。
その意味でRK4は、
高精度・汎用性・実装の簡潔さのバランスが非常に良い
という大きなメリットを持っています。
✨ まとめ
- 万有引力による運動方程式は2階常微分方程式である。
- RK4を適用する場合は、これを位置と速度による1階連立常微分方程式へ変換する。
- RK4は4次精度を持ち、滑らかな万有引力系を高精度に積分できる。
- RK4は非シンプレクティックなので、理論上、無限に長い積分でエネルギー誤差が完全に保存されるわけではない。
- しかし、初期値と時間刻みを適切に設定すれば、実用的な時間スケールでは非常に高い精度を維持できる。
- 実際に $e=0.95$、周期75年の軌道を50周、3750年相当積分したところ、軌道長半径の再現率は100.000%、離心率は100.002%だった。
- したがって、数値積分法を選択するときは「RK4かシンプレクティックか」という教条的な二択ではなく、対象とする物理系・時間スケール・時間刻み・必要精度を実測して判断することが重要である。
- スパコンを必要とするような極端な長期間・大規模多体系シミュレーションでなければ、一般的なPCで実用的な天体シミュレーションを行う場合、精度・実装の容易さ・計算量のバランスからRK4で十分なことが多い。
特殊なソルバーに頼り切る前に、まずは手元の数理を見つめ直し、初期値と時間刻みを極限まで研ぎ澄ましてみる。
それだけで、RK4というシンプルで信頼性の高い道具は、手元のPCの中に美しい小宇宙を立ち上げることができます。
作者
GitHub: https://github.com/NAS6mixfoolv
X（旧Twitter）: https://x.com/NAS6_oxo
作者HP: https://nas6.net
気に入っていただけたら GitHub に ⭐ をいただけると嬉しいです！