Einstein field equations (EFE) $N$-body simulator I made a $N$-body simulator and it works well, but it uses Newton's gravitational equation, which is nice, but I want it to simulate Einstein field equations. Speed of gravity should be simple enough, but how would I implement perihelion precession and the speed of an object increasing it's gravitational effect?
 A: There is no short answer to your question.
They are complicated since they are non linear and for that reason one has to use methods of approximation and last but not least they takes pages and pages of equations, proper surface integrals handling, and mastering general relativity.
This equations of n-body problem in GR are described in the following Einstein's papers :

*

*The Gravitational Equations and the Problem of Motion by A. Einstein, L. Infeld and B. Hoffmann(Annals of Mathematics, Second Series, Vol. 39, No. 1 (Jan., 1938), pp. 65-100)

*The Gravitational Equations and the Problem of Motion. II by A. Einstein and L. Infeld (Annals of Mathematics, Second Series, Vol. 41, No. 2 (Apr., 1940), pp. 455-464)

*On The Motion Of Particles In General Relativity Theory by A.Einstein (Canad. J. Math. 1, 209-241)

*On a Stationary System With Spherical Symmetry Consisting of Many Gravitating Masses Albert Einstein (Annals of Mathematics, Second Series, Vol. 40, No. 4 (Oct., 1939), pp. 922-936)

*and I would add this nice slide lecture The Two-Body Problem in General Relativity by Thibault Damour (Institut des Hautes Etudes Scientifiques (Bures-sur-Yvette, France))

As I started to give details here and here one must continue to the next steps by choosing the approximation method (described in the above papers), with the field of equation in mind and knowing which is which and what goes where.
I always advice people to read Einstein's book The meaning of Relativity in order to get into the correct spirit of General Relativity.
Only 36+10+32+15+...=93pages to read, but if we are satisfied with the first level of approximation which comes right after Newton's equations here's how it goes after doing all the math for 2 bodies:
$$ \ddot{\eta}^m - \overset{2}{m} \frac{\partial(1/r)}{\partial \eta^m} = \overset{2}{m} \left\{ \left[ \dot{\eta}^r \dot{\eta}^r  + \frac{3}{2}\dot{\xi^r} \dot{\xi^r} - 4\dot{\eta}^r\dot{\xi^r} - 4 \frac{\overset{2}{m}}{r} - 5 \frac{\overset{1}{m}}{r} \right] \frac{\partial}{\partial \eta^m} (1/r) + [4 \dot{\eta}^r(\dot{\xi}^m - \dot{\eta}^m) + 3\dot{\eta}^m \dot{\xi}^r -4\dot{\xi}^r\dot{\xi}^m ] \frac{\partial}{\partial \eta^r}(1/r) + \frac{1}{2} \frac{\partial^3 r}{\partial \eta^r \partial \eta^p \partial \eta^m}\dot{\xi}^m \dot{\xi}^p \right\}$$

The equations of motion for the other particle are obtained by replacing $\overset{1}{m}$, $\overset{2}{m}$, $\eta$, $\xi$ by $\overset{2}{m}$, $\overset{1}{m}$, $\xi$, $\eta$ respectively.

$\eta^r$ is the position of $\overset{1}{m}$, while $\xi^r$ is the postion $\overset{2}{m}$, r is the usual distance between the two objects. It's assumed the usual $c = G = \epsilon = \mu = 1$ (so you need to place back this constants to get proper units), the dots over the position vectors represent how many times do we derivate in respect to time (one dot = 1 time = velocity, 2dot = acceleration)
Note: here $m$ stands for $M$ where $M = \lambda^2 m $, and $\tau = \lambda x^0$, but we only noted $m$ for $M$ and $x^0$ for $\tau$. This is done to make the notation easier to understand; of course we still need to find $\lambda$ first.
