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.)