３Ｄプログラミング入門講座・その３１：物理演算その６
一般相対性理論：シュワルツシルト解の利用法


計量 ( ＋, －, －, － ) のシュワルツシルト解（球対称真空解）の線素は以下の通りです。

ds^2 = (1 - 2m/r) c^2 dt^2 - dr^2/(1 - 2m/r) - r^2 dθ^2 - r^2 sin^2θ dφ^2

ここで m = GM/c^2 （rs/2、rsはシュワルツシルト半径）です。

多くの解説では線素の導出に重点が置かれていますが、実際にシミュレーションでどう利用するかの情報は少ないため、
本稿では実装寄りの解説を行います。



シュワルツシルト解から1PN近似を導出する

赤道面（θ = π/2）を仮定します。 保存量は以下の2つです。

・エネルギー保存： (1 - 2m/r)(dt/ds) = l （定数）
・角運動量保存： r^2 (dφ/ds) = h （定数）

これらから導かれる u = 1/r の軌道方程式は

d^2u/dφ^2 + u = m/h^2 + 3m u^2

となります。右辺の 3m u^2 が相対論的補正項です。

補正項の物理的意味

加速度の形で書くと次のようになります。

a = ‐ (GM / r^2) × (1 + 3 (v/c)^2)

瞬間的な1PN補正因子は S = 3 (v/c)^2 です。

本ページでは、円軌道近似に基づく実装用の簡略化した1PN補正式を用いています。

2.5PN（重力放射減衰）の考慮

2.5PN(軌道平均近似) : 重力放射減衰 L = (-32/5) × (G^3 M m (M+m)) / (C^5 r^4) × v　[m/s^2]

シミュレーション用瞬間補正 : L' = L / P　（1周平均の減衰を各時刻へ経験的に配分）

※これは軌道平均式をシミュレーションへ適用するための実装上の近似です。

以下はニュートン万有引力に1PN,2.5PNを適用した加速度の計算例です。

実装例（C++風疑似コード）



double dt = deltaT;          // ΔT
double M = AttractorMass;    // 重力源質量
double m = SatelliteMass;    // 衛星質量
Vector3 r = AttractorPos - SatellitePos; // 中心方向ベクトル
double dist = r.length();                // 半径の大きさ
double v = velocity.length();            // 速度の大きさ
double v2 = v * v;
double r2 = dist * dist;

// ニュートン万有引力
double newton = G * M / r2;

double G3 = pow(G, 3);
double c2 = c * c;
double c5 = pow(c, 5);
double r4 = pow(dist, 4);

// 1PN補正
double pnFactor = 1.0 + 3.0 * (v2 / c2);

// 2.5PN（重力放射減衰）補正
double P = OrbitalPeriod * 86400 * 365.2425;  // 周期（秒）
double L25 = (-32.0/5.0) * (G3 * M * m * (M + m)) / (c5 * r4) * v;
double pn25Factor = L25 / P;

Vector3 accel_newton = r.normalized() * newton * pnFactor;   // 1PN:中心方向
Vector3 brake_pn25   = velocity.normalized() * pn25Factor;   // 2.5PN:速度逆方向への減衰

Vector3 acceleration = accel_newton + brake_pn25;

//###### RK4等の積分器側 ######

Vector3 velocity = velocity + acceleration * dt; // ここでdtがかかってくることに注意
この形で4次のルンゲクッタ法などに簡単に組み込むことができます。

注意点



本ページの目的は、シュワルツシルト解から得られる実用的な1PN,2.5PN補正式を
数値シミュレーションに組み込むための指針を示すことです。

一般相対論の理論を詳しく解説した資料は数多くありますが、
実際にシミュレーションへ組み込む手順をまとめた資料は意外と少ないように感じます。
本ページが、その橋渡しになれば幸いです。



S = × (v/c)^2
SA = -0.5 × (v/c)^2：特殊相対性理論
SB = 3.0 × (v/c)^2：一般相対性理論
SC = 2.5 × (v/c)^2：一般相対性理論＋特殊相対性理論
S0 = 0.0 × (v/c)^2：古典万有引力
つまり
SB = -6SA
SC = -5SA
ちなみに
SD = -1.0 / (v/c)^2 × (v/c)^2：等速直線運動：無重力
SE = -2.0 / (v/c)^2 × (v/c)^2：双曲線軌道：斥力反重力

F=-(GmM/r^2)(1+S)
S=(3.0)(v/c)^2これを例えば
S=(2.5e6)(v/c)^2として
大げさに相対論補正項を設定して
その振る舞いを見る



ちゃんとこういう軌道になっています

Source Code




















＜＜ prev　重力多体問題をマネージドクラス化する　：　グレゴリオ暦日付計算ライブラリ化　next＞＞




目次に戻る