计算科学 - 水星进动 Python 模拟 - 吾爱随笔录

我试图根据扰动的解决方案来模拟水星的进动：

\frac{1}{r} = \frac{m}{B^{2}} (1 + e \cos ϕ + 3 \frac{m^{2}}{B^{2}} (1 + e ϕ \sin ϕ + e^{2} (\frac{1}{2} - \frac{1}{6} \cos 2 ϕ))),

$\frac{1}{r}=\frac{m}{B^{2}}(1+e\cos\phi+3\frac{m^{2}}{B^{2}}(1+e\phi \sin\phi +e^{2}(\frac{1}{2}-\frac{1}{6}\cos2\phi))),$ 在哪里

e

$e$ 是汞的偏心率，B 是常数，它是单位质量的角动量。这是我在 Python 中的代码：http: //sprunge.us/TFhU

在代码中，我必须乘以一个大数（SF），以便稍微可见进动（或者你可以制作更大的 SF 值，那么进动会更容易看到，但它会扭曲轨道）。我想如果我用这个作为水星进动的插图是可以的，但是，我希望能够获得被称为每世纪 43 角秒的进动率的值，有没有人知道如何在代码或任何人的python代码做得更好？（我还通过了一个基于这两个方程的模型：

r_{n + 1} = 2 r_{n} - r_{n - 1} + Δ τ^{2} (- \frac{G M}{{r_{n}}^{2}} + \frac{B^{2}}{{r_{n}}^{3}} - \frac{3 G M B^{2}}{{r_{n}}^{4}}), ϕ_{n + 1} = ϕ_{n} + Δ τ \frac{B}{[\frac{1}{2} (r_{n + 1} + r_{n})]^{2}},

$r_{n+1} = 2r_{n}-r_{n-1}+\Delta \tau^{2}(-\frac{GM}{{r_{n}}^{2}}+\frac{B^{2}}{{r_{n}}^{3}}-\frac{3GMB^{2}}{{r_{n}}^{4}}),\phi_{n+1}=\phi_{n}+\Delta\tau\frac{B}{[\frac{1}{2}(r_{n+1}+r_{n})]^{2}},$ 其中 B 与上面相同，但我无法让它工作。）

（更新：我使用 Vpython 和 Leapfrog 方法来完成它，模拟结果很好。这是代码并根据需要修改它， http ://sprunge.us/ehgS （在 linux 中，球体没有留下痕迹可能是因为 linux 中的 vpython 非常过时）。但是，我仍然需要找到一种方法来输出数据中的进动率 43"。）