What's the exact gravitational force between spherically symmetric masses? Consider spherical symmetric$^1$ masses of radii $R_1$ and $R_2$, with spherical symmetric density distributions $\rho_1(r_1)$ and $\rho_2(r_2)$, and with a distance between the centers of the spheres $d$. What is the exact force between them? I know point masses are a good approximation, but I'm looking for an exact formula. This would be useful for a gravity-simulation toy software.
--
$^1$ Assume for simplicity the idealization where tidal or centrifugal forces do not deform the spherical symmetric, i.e., the various mass parts are held in place by infinitely strong and rigid bonds.
 A: The Gauss's Law answer is certainly right, but for a different take on things, you might want to have a look at the Principia itself. Long before Gauss came along, Newton proved the result you seek with geometry.
Basically, you consider how a test particle ("corpuscle") outside a spherical shell of mass is pulled by different sections of that shell, assuming the inverse-square law of gravity. This is Theorem XXXI in Book I. Theorem XXXIV extends this to a test particle outside a solid sphere of uniform density, essentially by integrating. Theorem XXXV generalizes again to the force between two disjoint spheres of uniform density. Finally, Theorem XXXVI tells us that in the general case you seek, the force is the same as though all the mass were concentrated at the spheres' centers, using the principle of superposition.
A: In the event you are interested in some methodology or a transparent way on how to solve this problem, here is an analytic approach.
Let the two spheres be $S_1$ and $S_2$ radius $R_1$ and $R_2$ at centres C1 and C2 respectively. The distance $C1C2 =d\ge R_1+R_2$.  Work in spherical coordinates with the frame at C1:
1)Consider two point P1$(  r_1\cos\phi_1\cos\theta_1, r_1\cos\phi_1\sin\theta_1, r_1\sin\phi_1)$ 
and P2$( d+ r_2\cos\phi_2\cos\theta_2, r_2\cos\phi_2\sin\theta_2, r_2\sin\phi_2)$ in the volume of $S1$ and $S2$ respectively. 
2)  Now write the square of the distance P1P2
$r_{12}^2=(d+ r_2\cos\phi_2\cos\theta_2- r_1\cos\phi_1\cos\theta_1)^2+(r_2\cos\phi_2\sin\theta_2-r_1\cos\phi_1\sin\theta_1)^2+ (r_1\sin\phi_1- r_2\sin\phi_2)^2$
3)  Write the differential volume elements  for the two spheres at the points P1 and P2 in spherical coordinates:
$dV_1=\sqrt {g_1}dr_1d\phi_1d\theta_1$ and $dV_2=\sqrt {g_2}dr_2d\phi_2d\theta_2$
where $g_1$ and $g_2$ are the Jacobian determinants for spherical coordinates (not difficult to find, or look it up in a vector analysis book, or use $dV=r^2\cos\phi drd\phi d\theta$). Hence the differential masses of these volume elements are
$dm_1=\rho_1(r_1)dV_1$ and $dm_2=\rho_2(r_2)dV_2$
Now you can put all these together and write the  total magnitude of the force
$F=G\int_{V_1}\int_{V_2}\frac {\rho_1(r_1)\rho_2(r_2)dV_1dV_2}{r_{12}^2}$ 
and you need to do this integral with the transformations given above.
This is a very general formula giving the total force in the problem. One can see that if $d>>R_1+R_2$ this equation reduces to the gravitational attraction between two point masses at large distance $d$. One could end up with simpler integrals by using the symmetry in the line C1C2, and considering discs instead the general volume elements we have done in the above method
A: For spherically symmetric bodies the point-mass approximation is actually exact as a description of the centre-of-mass force on the spheres. This is due to Gauss's law in its gravitational form: the flux of the gravitational field through a concentric spherical surface is proportional to the mass inside it, or mathematically
$$\oint\mathbf{g}\cdot d\mathbf{S}=4\pi G M_\text{int}.$$
This can be seen to be equivalent to Newton's law of gravitation and it yields the point-mass formula when applied to a spherical mass.
A: If you are looking only for the Newtonian gravitational force, the other answers to this question are correct.  The spherically symmetric mass distributions can be replace by the total mass at the center of mass and then the Newtonian gravitational force can be computed for these two masses.
However, if you want EXACT calculations, you must use general relativity.  In general relativity it is also true that a spherically symmetric mass distribution can be replaced with the Schwarzchild black hole of the same mass at the same center.  However this is for a single spherical distribution.  It is probably not true for two spherical distributions or for two Schwarzchild black holes.  There is no know exact general relativity solution to the two black hole problem.  All we have is very good numerical simulations of general relativity for two black holes.
