general relativity - Modifying Newtonian gravity to fit observed precession of Mercury's orbit

The idea is to modify Newtonian gravity so that it fits measurements of orbits around the sun. For example the precession of Mercury's orbit unlike Newtonian $n$-body simulations.

I'm currently not using this in any serious simulation, but I'm just wondering because it'd simplify approximations a lot.

The modified formula for the force I found about half a year ago after googling was similar to this one: $$F =\frac{Gm_1m_2}{r^2} + \frac{Gm_1m_2B^2}{r^4} = \frac{Gm_1m_2}{r^2} (1+\frac{B^2}{r^2})$$

Similar means that I lost the link which I was unable to relocate and this is the only thing I wrote down elsewhere. $B$ stands for the dot product of velocity and unit vector pointing at the other object if I remember right.

This all sounds rather vague and arbitrary. Thus I'm asking whether anyone knows of something like this and the explanation behind it.

Answer

Here are couple of references that describe professional uses of a post-Newtonian formalism to model the planets and the Earth's Moon:

Standish, et al. "Orbital ephemerides of the Sun, Moon, and planets," Explanatory Supplement to the Astronomical Almanac (1992): 279-323.
The relevant equation is 8-1 on page 3.

Petit and Luzum (eds.), "IERS Technical Note No. 36, IERS Conventions (2010)," International Earth Rotation and Reference Systems Service, Frankfurt, Germany (2010). You'll want to look at chapter 10. This describes the Moon, but not the planets. The IERS's focus is the Earth.

Note that the IERS document does describes relativistic time scales. Mixing and matching the non-relativistic time scales we use to keep time on the surface of the Earth with a post-Newtonian formalism doesn't buy much. JPL (the first reference) uses its own relativistic time scale, T_eph. This is now very close to one of the officially defined time scales, Barycentric Dynamical Time (TDB). The IERS document uses a different relativistic time scale that doesn't tick with Earth-based clocks.

(Note: The acronym doesn't match the name. That's intentional. The official acronym is short for the French name "Temps Dynamique Barycentrique", regardless of what language you use for the full name.)

---

Edit: Regarding the simplified version specified in the question

My old, old copy of Marion, Classical Dynamics has something very close to the formula in the question. It extends the Newtonian relation

$$ \frac{d^2}{d\theta^2}\left(\frac 1 r\right) + \frac 1 r = -\frac m {l^2}r^2 F(r) = GM \frac {m^2}{l^2}$$

to

$$ \frac{d^2}{d\theta^2}\left(\frac 1 r\right) + \frac 1 r = GM \frac {m^2}{l^2} + \frac {3GM} {r^2 c^2} $$

from which one the force modified $F(r)$ can be expressed as

$$F(r) = -\frac {GMm}{r^2}\left( 1 + \frac{3 \, l^2}{m^2c^2 r^2} \right)$$

Using $\vec l \equiv \vec r \times (m \vec v)$, this can be rewritten as

$$F(r) = -\frac {GMm}{r^2} \left( 1 + \frac{3 \, ||\vec r \times \vec v||^2}{c^2 r^2} \right) $$