What is the exact gravitational force between two masses including relativistic effects? Is there a closed-form formula for the force between two masses $m_1$ and $m_2$ if relativistic effects are included? My understanding is that the classic formula $G \frac{m_1 m_2}{r^2}$ is just an approximation (which is probably good enough even for going to the Moon), but what would the correct formula be according to the general theory of relativity? Does a closed formula even exist? For example for an idealized situation of two uniform spherical masses?
 A: Yes the non-relativistic equation of motion 
$$\frac{dp}{dt} = F = -\frac{G m_1 m_2}{r^2}$$
for a single particle interacting gravitationally with other is only valid for velocities much smaller than the speed of light and for small enough masses. Notice the minus sign in the expression for the non-relativistic force F --gravitation is attractive--.
First, general relativity is a (geo)metric theory. There is no gravitational forces in general relativity. In general relativity, bodies affected only by gravitation are moving freely, but in a curved spacetime

(source: twimg.com) 
The general relativistic equation of motion for a single body is the geodesic equation
$$\frac{DP^\mu}{D\tau} = 0$$
Note the zero at the right, which is a consequence of the absence of gravitational forces in general relativity. The greek indices run over 0,1,2,3 --the spacetime coordinates-- and summation convention is being used. $P^\mu$ is the four-momentum, $\tau$ the proper time and $D$ denotes the covariant derivative, which includes the effects due to spacetime curvature
$$\frac{DP^\mu}{D\tau}
=
\frac{dP^\mu}{d\tau}
+
\Gamma_{\nu\lambda}^\mu U^\nu P^\lambda ,$$
where $\Gamma_{\nu\lambda}^\mu$ are the Christoffel symbols and $U^\nu$ the four-velocity.
Second, the field theory of gravity provides a non-geometrical description of gravity. There are gravitational forces in the field theory. In this theory, bodies affected only by gravitation are moving in a flat spacetime --sometimes this is named the flat spacetime approach to gravitation-- but feel a gravitational force associated to gravitons

(source: twimg.com) 
The field-theoretic equation of motion for a single body in a gravitational field is the Kalman equation
$$A_\mu^\nu \frac{dP^\mu}{d\tau} = F^\nu  = - B_{\mu\lambda}^\nu U^\mu P^\lambda$$
where $F^\nu$ is the exact gravitational force with
$$A_\mu^\nu = \left( 1 - \frac{1}{c^2} \psi_{\lambda\gamma} U^\lambda U^\gamma \right) \eta_\mu^\nu - \frac{2}{c^2} \psi_{\mu\gamma} U^\gamma U^\nu + \frac{2}{c^2} \psi_\mu^\nu$$
and
$$B_{\mu\lambda}^\nu = \frac{2}{c^2} \psi_{\mu,\lambda}^\nu - \frac{1}{c^2} \psi_{\mu\lambda}^{,\nu} - \frac{1}{c^2} \psi_{\mu\lambda,\gamma} U^\gamma U^\nu$$
Here $\psi_{\alpha\beta}$ is the gravitational field potential and the comma denotes the ordinary flat spacetime partial derivative.
A: As in Electrodynamics, the forces become retarded, the corresponding equations become complicated and include radiation too.
A: The expression that is actually used by JPL to calculate/approximate the relativistic effects under given promises comes from the "post-Newtonian expansion" at the 1PN-level and, including the classical Newtonian acceleration, it looks like:
$$\frac{d\bar{v}}{dt}=-\frac{GM}{r^2}\left(1-\frac{4GM}{rc^2}+\frac{v^2}{c^2}\right)\hat{r} +\frac{4GM}{r^2}\left(\hat{r}\cdot \hat{v}\right)\frac{v^2}{c^2}\hat{v}$$
So they are basically introducing one component that is "repulsive inverse cube gravity" and two velocity dependent terms. This expression is also availabe at the 3PN-level, see this post. Note that this is only an approximation that is able to reproduce the "anomalous perihelion shift" in the weak field limit. The extra repulsive inverse cube term causes very strange behaviour if you try to use it in the strong field limit. Here are a few simulations I made where the green circle denotes the Schwarzschild radius and the red circle denotes the radial distance to the "innermost stable circular orbit." It is easy to see the bouncing caused by the repulsive inverse cube term:

Note: This is the expression for the case of infinite mass ratio. If you cannot use the infinite mass ratio approximation, the expression would be more complicated I think.
A: Worse than electrodynamics, general relativity is non-linear, in the sense that the field from multiple sources is not just the sum of fields from each isolated source. Even the simple case you are asking about, which is the two-body problem in general relativity, has not been solved exactly.
An even simpler case is the limit $m_2 \to 0$. In that case only $m_1$ affects the geometry of spacetime, and $m_2$ follows a geodesic in that spacetime. This has been solved exactly. The linked article gives details.
[Addendum] To directly answer the original question for the limit $m_2 \to 0$:
Of course if $m_2=0$, then the force between the two masses is $0$. But the thing about classical gravity is that the acceleration due to gravity is independent of mass.  (This is the equivalence principle and is actually one of the starting points of the theory of general relativity.) So it still makes sense to ask what would be the acceleration of a negligible mass (aka, a test body) due to the gravity of another mass. The classical answer is $a_2 = \frac{Gm_1}{r^2}$.
In general relativity, the acceleration of a test body due to the gravity of a single spherical, homogenous, non-rotating mass is given exactly by the Schwarzschild solution, which link you can consult for details. The result of which is that
$$\ddot{r} = -\frac{Gm_1}{r^2} + r\dot{\theta}^2 - \frac{3Gm_1}{c^2}\dot{\theta}^2$$
$$\ddot{\theta} = -\frac{2}{r}\dot{r}\dot{\theta}$$
[CORRECTED AND SIMPLIFIED Jan 2]
where $r$ and $\theta$ are polar co-ordinates centred on the gravitating mass, the dots represent differentiation by the proper time of the test body.
So the first term is just the classical radial acceleration $-\frac{Gm_1}{r^2}$. The terms $r\dot{\theta}^2$ and $-\frac{2}{r}\dot{r}\dot{\theta}$ are the classical centrifugal and Coriolis acceleration for polar co-ordinates.
What's not classical is the extra term $\frac{3Gm_1}{c^2}\dot{\theta}^2$. Finally, there is the fact that differentiation is with respect to proper time of the test body. Different test bodies will experience time differently. They can be related by:
$$(1 - \frac{r_s}{r})\dot{t}^2 - \frac{\dot{r}^2}{(1 - \frac{r_s}{r}) c^2} - \frac{r^2\dot{\theta}^2}{c^2} = 1$$
The constant $r_s = \frac{2Gm_1}{c^2}$ is introduced for simplicity. It is called the Schwarzschild radius of the gravitating mass. 
Here the co-ordinate $t$ is introduced as a reference time, so $\dot{t}$ is the rate of change of reference time with respect to proper time of the test body. For a distant ($r \to \infty$), stationary ($\dot{r}=0, \dot{\theta}=0$) test body, this becomes $\dot{t} = 1$, so the reference time can be interpreted as the time measured on a distant, stationary clock.
In the classical case of course every body experiences the same time, but it can also be compared to the special relativistic case, where the equation would be:
$$\dot{t}^2 - \frac{\dot{r}^2}{c^2} - \frac{r^2\dot{\theta}^2}{c^2} = 1$$
So what's new in general relativity is the factor $(1 - \frac{r_s}{r})$. To get some idea of scale, for the Earth, $\frac{r_s}{r}$ is about one and a half parts per billion on the surface of the Earth. (Note the Schwarzschild solution is only applicable outside the gravitating body.)
This is exact only for $m_2 \to 0$, but it remains a very good approximation as long as $m_2$ is much smaller than $m_1$, such as for a planet orbiting a star.
A: I think if we consider a much bigger mass then the force becomes relativistically: 
$$ F=-\frac{GMm}{r^2}+\frac{4G^2M^2m}{r^3c^2}$$
