# Precession of Mercury (Python simulation)

I was trying to simulate the precession of Mercury based on the perturbed solution, and my questions about its implementation in python can be seen here: https://scicomp.stackexchange.com/questions/19778/precession-of-mercury-python-simulation.

And as I was searching on-line, I saw many papers were using post-Newtonian method to do it, but the formula involved in PN always have very complicated expression and I didn't really find any code implementation, why don't they just use the perturbed solution(the solution in the link above) coming from Lagrangian(the geodesic of time-like particle) instead of using those complicated expression? Can someone explain?

(EDIT: here is a blog that talks about the post-Newtonian method approach: https://astrokode.wordpress.com/2014/05/03/the-precession-of-mercurys-perihelion-simulation/)

(Update: Eventually I used Leapfrog method and adding some correction to the Newtonian gravity to do the simulation in Vpython, the code can be found in the first link. And I assume the reason why people use PN method is simply because in real life you don't get static and nice spherically symmetric objects, and just use the perturbed solution in schwarzschild metric would not be accurate in a long time simulation.)

• Might Computational Science be better suited for the programming aspect of this question? Commented May 29, 2015 at 1:02
• @KyleKanos I didn't know that exists, just posted one there. Thanks for the info
– Sam
Commented May 29, 2015 at 1:15
• Our comp-phys is somewhat limited in scope (not about implementations). Typically, SE does not allow for cross-posting; your options would be either deleting the question & re-asking it on Computational Science or flagging it for moderator attention. I personally think that the last question (about PN corrections) should be on-topic here and can be left alone (it really stands alone as a separate question from the first bit anyways), so maybe you'd want to remove the first bit & post that as a separate question on Computational Science? Commented May 29, 2015 at 1:15
• @KyleKanos hope it's better now.
– Sam
Commented May 29, 2015 at 1:25
• Mercury's relativistic precession is 43 arc seconds per century. That is tiny. I truly doubt you'll see such a subtle effect using leapfrog. The blog you cited uses IAS15, a 15th order integrator specifically designed for gravitational dynamics. Commented Jun 6, 2015 at 11:51

It's quite easy to get the necessary precession for Mercury if you take gravitational and kinematic time dilation / length contraction into account.

The Newtonian acceleration vector is: $$$$\vec{g}_n = \frac{\hat{d} G M}{{\lvert\lvert \vec{d} \rvert\rvert}^2}.$$$$

One important value is closely related to the kinematic time dilation: $$$$\label{eq_kinematic} \alpha = 2 - \sqrt{1 - \frac{\lvert\lvert \vec{v}_{o}\rvert\rvert^2}{c^2}}.$$$$ Another important value is the gravitational time dilation: $$$$\beta = \sqrt{1 - \frac{R_{s}}{\lvert \lvert \vec{d} \rvert \rvert}}.$$$$

Finally, the semi-implicit Euler integration is: \begin{align} \vec{v}_{o}(t + \delta_t) &= \vec{v}_{o}(t) + \delta_{t} \alpha \vec{g}_n, \\ \ell_{o}(t + \delta_t) &= \ell_{o}(t) + \delta_{t} \beta \vec{v}_{o}(t + \delta_t). \end{align}

I won't pass on my code quite yet, but here's a beginner's C++ source code that calculates strictly Newtonian orbit. Adding in the required code isn't all that tough, but I just don't know if my code has bugs, so I need independent validation:

https://github.com/sjhalayka/mercury_orbit_glut

Let me know if you have any questions!

... and yes, I have already submitted this to the journal.