JavaScript で一般相対論：シュワルツシルト解の 1PN/2.5PN を万有引力に掛けて軌道を描く

Last updated at 2026-08-11Posted at 2026-08-11

JavaScript で一般相対論：シュワルツシルト解の 1PN / 2.5PN を万有引力に適用する

一般相対性理論の軌道は難しいと思われがちですが、
実は ニュートン万有引力に「(1+S)」を掛けるだけで、
かなり直感的に相対論的軌道を再現できます。

この記事では、以下の HTML ページで実装している

F = -(GM/r²)(1+S)
S = α (v/c)²

という補正モデルを使い、
1PN（一般相対論）・2.5PN（重力波放射減衰）を JavaScript で可視化する方法を紹介します。

📐 シュワルツシルト解から 1PN を取り出す

赤道面（θ = π/2）を仮定すると、
シュワルツシルト解の軌道方程式は次の形になります。

$$
\frac{d^2u}{d\phi^2} + u = \frac{m}{h^2} + 3mu^2
$$

右辺の 3mu² が一般相対論の補正項です。

これを加速度の形にしようと試みます：

ここでは弱い重力場・低速近似（v≪c）の下で
ds≈c dt と置き、接線速度を軌道速度 ( V ) とみなして
整理しています。

この式は
$$\frac{d^2u}{d\phi^2} + u = \frac{m}{h^2}$$
が万有引力の式に対応しイメージしてみて
r半径,G万有引力定数,M質量,V軌道速度,c光速度,S相対論補正項
$$a=-(\frac{GM}{r^2})(1+S)$$
の形であるから、ここで
$$h=\frac{r^2d\phi}{ds}$$,$$m=\frac{kM}{c^2}$$,$$u=\frac{1}{r}$$,
であって
相対論補正項$$3mu^2$$を検証すると、
$$\frac{m}{h^2}$$でくくり
$$u''+u=(\frac{m}{h^2})(1+\frac{3mu^2}{\frac{m}{h^2}})=\frac{m}{h^2}(1+S)$$の形にして、Sを調べます
$$\frac{3mu^2}{\frac{m}{h^2}}=3u^2h^2=3(\frac{1}{r^2})(\frac{r^2d\phi}{ds})^2$$
$$=3r^2((\frac{dt}{ds})(\frac{d\phi}{dt}))^2=3r^2(\frac{d\phi}{dt})^2(\frac{1}{c^2})$$
$$=3(\frac{rd\phi}{dt})^2(\frac{1}{c^2})$$
$$(\frac{rd\phi}{dt})$$は円の接線方向の速度だからそれを軌道速度Vとみれば$$3(\frac{V}{c})^2$$と解けて
つまり相対論補正項$$S=3(\frac{V}{c})^2$$だから
ユークリッド幾何で
シュワルツシルト時空の重力による加速度の式は
r半径,G万有引力定数,M質量,V軌道速度,c光速度,S相対論補正項
$$a=-(\frac{GM}{r^2})(1+S)$$
$$=-\frac{GM}{r^2}\left(1 + 3\left(\frac{V}{c}\right)^2\right)
$$
と解ける
この 1PN 補正項は、ニュートン項（∝ 1/r²）に対して
おおむね 1/r³ のオーダーの大きさです。
一方、次の 2PN 以降の補正は 1/r⁴ 以下のオーダーになるため、
多くの実用的な場面（水星の近日点移動など）では
1PN だけを考慮すれば十分であることが多いです。
ただし、強い重力場（ブラックホール近傍など）や、
極めて高い精度が要求される場合は、より高次の項も重要になります。つまり r が極端に小さい場合などです

📉 2.5PN（重力波放射減衰）

連星が重力波を放出すると軌道エネルギーが減少し、
軌道は徐々に縮んでいきます。

軌道平均近似では次の式になります。

$$
L = -\frac{32}{5}\frac{G^3 M m (M+m)}{c^5 r^4} v
$$

シミュレーションではこれを 1周あたりの減衰量を各時刻に配分する形で使います。つまりLを軌道周期で割る必要があります。
そうして計算された値を軌道速度の逆方向の成分として加算(つまりブレーキ)します。

🧪 JavaScript 実装（抜粋）

以下は実際の HTML ページで使っている実装の一部です。

// 1PN 補正
N6LRngKt.prototype.ToSchwartz = function (v) {
    return gsval * v * v; // S = α(v/c)^2
};

// 加速度計算
N6LRngKt.prototype.GetA = function (r, m, mr, v) {
    var g = this.CNST_G * (m / mr / mr);
    g = g * (1.0 + this.ToSchwartz(v)); // 1PN
    return g / this.CNST_C;
};