Why is the Force of Gravitational Attraction between two “Extended” bodies proportional to the product of their masses? Newton’s Law of gravitation states that force of attraction between two point masses is proportional to the product of the masses and inversely proportional to the square of the distance between them. I know that the force of attraction between two spheres turns out to be of the same mathematical form as a consequence of Newton’s law. But I am not able to prove how the force between any two rigid masses is only proportional to the product of their masses (as my teacher says) and the rest depends upon the spatial distribution of the mass. So $F$ is ONLY proportional to $Mmf(r)$ where $f(r)$ maybe be some function based on the specifics of the situation.
 A: How does this sound?
Let's say you have two extended bodies, A and B, each made up of a number of particles. Let's consider the force on a particle in A, call it P, due to body B. Each particle in body B exerts a force on P that is proportional to that B particle's mass. The sum of such forces gives the net force on P due to B. Now let's suppose we double the mass of each B particle. That would result in doubling the force on P. But doubling the mass of each B particle is doubling the total mass of body B. So doubling the total mass of body B doubles the force on P. Hence, the force due to body B is proportional to B's total mass.
To get the total force due to body B on the body A, we have to add up the forces on all of A's particles. Again, if we doubled each particle's mass, the total force on A would double, as would its total mass. Hence that total force is proportional to A's total mass.
A: As is stated in the previous answer huge non-spherical objects are rarely (never) seen in Nature so let's consider this as a purely theoretical enterprise. Let's go!
Consider two bodies of mass (let's assume they are rigid so they can't be deformed and tidal effects omitted). The gravitational force vector a massive body exerts on another body originates in its center of gravity (CG). Also, the point where it gets a grip on the other body is the other body's CG. For most celestial bodies (which are spherical symmetric) as well as for a number of other bodies where symmetry is involved (think of an ellipsoid), the center of mass (CM) and CG coincide, but in general, this is not the case. In fact, there is for these bodies not one CG, but they lie on a line piece on which the CM is situated. Where the CG is situated depends on where the other body is to be found. For example on a cube Earth: 

The direction of the gravity you feel when walking on it doesn't point most of the time towards the CM. In the article is written:

“… Gravity on the surface wouldn’t generally point toward the exact center of the [cube] Earth anymore.”

And the center of this cube Earth is where the CM is situated. 
The same holds for a barbel formed mass. Although some symmetry is involved the CG lies on a line piece on the principal axis of rotation with the smallest moment of inertia (in the barbel's length) going through the CM. The endpoints of the line piece lie a bit before the center of the two balls, as seen from the CM, and on the same distance from it (unless the balls are different in size, in which case the CM shifts from the middle and doesn't lay in the middle of the line piece anymore).
And also in this case the position of the CG depends on the position of the other body.
Now, what does this all mean? It means that if two rotating bodies of mass (with masses $M$ and $m$ gravitate towards each other in a straight line, the forces of gravity lay on the line connecting the two CGs that lay on one of the endpoints of the line of CGs (dependent on their relative orientations). This means they make a torque (a pseudovector produced by a cross product) come into existence on both bodies: 
$$\vec{{\tau}_{1net}}=\vec{F_{gnet}}\times \frac{1}{2}(\overrightarrow{CM_1}-\overrightarrow{CG_{1max}})=\lVert{\vec{F_{gnet}}}\rVert \lVert{\frac{1}{2}(\overrightarrow{CM_1}-\overrightarrow{CG_{1max}}})\rVert \sin{{\theta}_1}$$
$$\vec{{\tau}_{2net}}=\vec{F_{gnet}}\times \frac{1}{2}(\overrightarrow{CM_2}-\overrightarrow{CG_{2max}})=\lVert{\vec{F_{gnet}}}\rVert \lVert{\frac{1}{2}(\overrightarrow{CM_2}-\overrightarrow{CG_{2max}}})\rVert \sin{{\theta}_2}$$
The factor $\frac{1}{2}$ appears in front of the vectors because their magnitudes vary between zero (when the CMs and CGs coincide) and the maximum values (when the CMs and CGs are the most far apart). This is the case on both sides of the CM, but on one side the $\vec{F_g} s$ have bigger magnitudes, which is why I write $\vec{\tau_{net}}$ and $F_{gnet}$, which I will call just $\vec{\tau}$ and $\vec{F_g}$ in what follows.
When the bodies have an initial minimum angular rotation the bodies make a whole rotation, while they rotate back and forth when they rotate below this minimum angular rotation, and no rotation at all results when the initial angular momentum is zero and the line piece between the two CGs is perpendicular to the line between the two CMs of the bodies.
The angles ${\theta}_1$ and ${\theta}_2$ are the angles between the concerning vectors. Their maximum values increase when the bodies get closer. The torque vectors are perpendicular to the plane containing $\vec{F_g}$ and $(\overrightarrow{CM_1}-\overrightarrow{CG_{1max}})$ or $(\overrightarrow{CM_2}-\overrightarrow{CG_{2max}})$ (if these are parallel, no torques are present, because in that case $\sin{\theta}_1$ and $\sin{\theta}_2$ are zero) and rotate around the axis connecting ${CM}_1$ and ${CM}_2$.
Now a torque makes a body rotate in the plane perpendicular (or rotating back and forth as we saw) to the plane just mentioned and are also given by:
$$\vec{{\tau}}_M=I_M\frac{\vec{{d\omega}}_M}{dt}$$
$$\vec{{\tau}}_m=I_m\frac{\vec{{d\omega}}_m}{dt},$$
where $I_M$ is the moment of inertia of the body of mass $M$ and  $I_m$ the moment of inertia of the body of mass $m$ (different momenta of inertia $I$, depending on the form of the mass, can be calculated or looked up) and $\frac{\vec{d\omega}}{dt}$ the time derivative of the angular velocity. It might be clear that the time derivative of the angular velocity (pseudo)vector ($\vec{\omega}=\vec{v}\times \vec{l}$) lies on the same line as the torque vector because $I$ is a scalar, i.e. a positive number. Let's assume the torques makes the bodies rotate only around the principal axes with the highest moment of inertia (the rotation can be around any axis, but the principle is the same).
When the bodies initially don't rotate, and the CMs and the CGs of both bodies coincide or all lie on one line (which is the case if the lines of CGS are parallel or orthogonal, so no torque is present) then the bodies just accelerate linearly to each other with a force $F=G\frac{Mm}{r^2}$, where $r$ is the distance between both CMs (or CGs).
If this is not the case the bodies make each other rotate. After each full rotation of each body, they have the same angular velocity.  
So we have the mutual linear acceleration caused by the force $F=G\frac{Mm}{r^2}$, in which $F$ is the force component of the force pulling on the CMs and $r$ the distance between the two CMs. This linear acceleration is minimal when the total rotational acceleration is maximal, and vice-versa.  
This linear acceleration is attenuated periodically by the variable (but periodic) angular rotation of both bodies. The linear acceleration grows, diminishes, grows, diminishes, etc. (because the rotational energy of the bodies varies periodically).
When the bodies are very far apart there will (approximately) only be linear acceleration, because the torques go to zero. But the linear acceleration also goes to zero when they are far apart so $f(r)$ reduces to one, so you can reduce the force formula to $G\frac{Mm}{r^2}$, with $M$ and $m$ considered as point particles. The torques and linear acceleration will not have the same ratio when the distance between the bodies increases.
The torques of the bodies (making their rotation vary between a maximum and minimum value) are just as the force $F_g$ (producing the linear acceleration) dependent on a squared distance ($F\propto{\frac{1}{{r_F}^2}}$ and $\tau\propto{\frac{r_{\tau}^2}{r_F}^2}$) so initially, the linear acceleration is very small, just as the torques. The angles involved in the torques go to zero though when the bodies approach infinity and so does $\sin{\theta}$ for both bodies.
So the attenuating function $f(r)\rightarrow 1$ when $r\rightarrow \infty$ so $F_g$ approaches the $G\frac{Mm}{r^2}$ form. When the bodies approach each other the ratio of the torque of the bodies and linear accelerating force is not equal at every distance between the bodies (see the previous alinea). So $f(r)$ is a periodic function (depending on the initial rotations of the bodies, the momenta of inertia of them both, and the varying torques) and it gives the linear acceleration a periodic component. This periodic variation is small when the linear acceleration is small (when they are far apart) and gets bigger when the linear acceleration increases (when they are getting closer). But because the ratio of the torque and F grows when the distance between $CM_1$ and $CM_2$ diminishes (e.g. when the distance gets half as big, the force becomes $\frac{1}{4}$ big, while the torque becomes more than $\frac{1}{4}$ as big, because in the two cross products definitions I gave above, a factor $\sin{\theta}$ is involved, which grows when the distance of the bodies gets smaller), the period of rotation in time get. Nevertheless, $f(r)$ still attenuates the linear acceleration periodically.
Given the necessary information, $f(r)$ can be calculated. Of course, we also have to include the stretching of the bodies because they are not truly rigid. This stretching is due to the rotations and the tidal forces (the last grows when the distance gets smaller), which can be calculated too. When the bodies stretch potential energy is given to the bodies and this diminishes the linear acceleration, but this effect I neglected (though it makes a small contribution to $f(r)$ and also approaches one when the distance approaches infinity) which is why I assumed them to be rigid.  
In most cases, there is no constant periodicity (i.e. after a certain number of rotations the initial relative position of the bodies arises again) which is the case when the ratio of the momenta of inertia is a non-rational (real) number, but nevertheless a periodicity is present.
Pffff....I think this is more than enough.
Just one more thing. I just realized that the CGs don't have to lay on a straight line piece (this is only the case when symmetry is involved), but in general, they lay on a curved line piece. As a first approximation, it will do though (like a first approximation in a multipole expansion).
A: The statement

...the force between Any two rigid masses is only proportional to the
  product of their masses

is not true in general, or at least it is misleading. The shapes of the mass distributions and their relative positions matter when computing the gravitational force.
It is true that once you hold constant the shapes of the mass distributions and their relative positions, then the force will be proportional to the product of the total masses of the bodies.
There are certain situations where treating two extended massive bodies as point sources can be exactly correct (in the context of Newtonian gravity). For a spherically symmetric mass distribution, the gravitational potential outside of it is the same as that arising from a point source of the same mass. This is an application of Gauss' law.
In general, one can build up an increasingly good approximation of the gravitational potential arising from a given mass distribution via a multi-pole expansion. . The leading-order term, which drops off least rapidly with distance (force $\propto r^{-2}$), is that of a monopole like what arises for a point mass or outside a spherically symmetric system. But a general mass distribution will have contributions from higher-order terms (dipole, quadropole, octopole...), all of which drop off increasingly rapidly with distance. As one considers two bodies at increasing separation, reducing them both to their monopole terms becomes increasingly more accurate.
Finally, the fact that the gravitational force of attraction on an extended body due to another body can vary with position is essential when considering phenomena such as tidal forces.
A: This is not true in the general case if the two masses are close to each other. That is, if you have 2 different objects of the same mass but with different shapes, the gravitational attraction between these 2 objects and a third mass that is close enough to the objects will depend on the specifics of the situation. For example, imagine a large dumbbell. In the middle between the two big masses at the end, the gravitational force of these masses is equal and opposite, so with a zero resultant. Only the small bar in between the two "counts". Compare that to a large plate of the exact same mass. In this case, all the mass "counts". Of course, if the objects have the same shape, but one has a mass double the other, then the forces will be doubled. 
If the objects between which you calculate attraction are far enough from each other however, this is a good approximation, irrespective of shape. To get the attraction you integrate the attraction between each small element of object 1 and each element of object 2. If you are not at the calculus level, you can say that you divide both objects in small pieces and add the contribution of attraction between all the pieces with each other. In that case, the 1/r^2 factor is approximatively a constant for all "force couples" and the inverse square law is a good approximation.
To address your comment in the question, changing the mass distribution but keeping shape the same can be equivalent to changing the distance, so this can change attraction. For example, in my dumbbell, make one of the dumbbell much more massive than the other one, then the attraction at the center changes direction.
As a last note, even if you use the center of mass of the objects to get the distance "r", changing mass distribution will affect attraction. The center of mass moves as a linear function of distance in the mass distribution, but the attraction moves as a function of 1/r^2. These cannot compensate each other exactly all the time. In my dumbbell example, you are at the center of mass in the middle between the two balls, yet attraction is zero. Now, transfer half of the mass of one sphere to the other. The center of mass moves towards the mass that gets heavier. Go to the new center of mass. The attraction is not zero, as you are now close to a larger mass and far from a small mass.
A: It is not true that the gravitational force between extended mass distributions depends only on the product of the total masses. It is true that the time averaged total force integrated over each body is $$\vec F = Gm_1m_2 \frac{\vec r_{12}}{r_{12}^3} ~.$$ However, unless both mass distributions are spherical, the attraction has higher moments. These higher moment forces cause the bodies to be stressed and to nonuniformly rotate or wiggle. Only for certain relative orientations these higher moment forces exactly cancel the internal stresses.
An example is the Earth-Moon system. The moon is deformed but it is almost at rest in the corotating frame. It only wiggles a little. Weirder is the rotation of Mercury. It has a slight permanent dipole deformation causing it to rotate in a tidal 3:2 resonance. See https://en.wikipedia.org/wiki/Mercury_(planet)#Spin-orbit_resonance.
A: If to be summarized in short - you need to apply and solve second Newton law equation for two-body problem :
$$ \vec F_G = \mu \, \vec r^{\,\prime \prime} $$
Where $\mu$ is two-body system reduced mass :
$$ 
\mu ={\cfrac {m_{1}m_{2}}{m_{1}+m_{2}}}
$$
Btw, it's interesting to note that reduced mass has reciprocal additive property :
$$ {\frac {1}{\mu }}={\frac {1}{m_{1}}}+{\frac {1}{m_{2}}} $$ 
Reduced mass helps to analize two-body problem as it were just 1 single body. And $\vec r$ is displacement between bodies.
That's why gravity force is proportional to mass product of both bodies. (I.e. product increases faster than sum of masses). Another way which is helpful to physical intuition is to check the moment of inertia of binary system :

Which is :
$$ 
I={\frac {m_{1}m_{2}}{m_{1}\!+\!m_{2}}}x^{2}=\mu x^{2}
$$
A: It is not true in general that the gravitational force of attraction between extended bodies is proportional to their masses. It happens that we usually deal with gravitational attraction between celestial bodies, and that celestial bodies above a certain size are almost invariably close to spherical (in consequence of the self gravity of the body). In the particular case of spherical bodies, the result is true as a consequence of Newton's shell theorem.
In the general case, simply note that the inverse square law of gravity is basically the same (up to the sign of charge) as the Coulomb law of electrostatics, and apply the argument of any number of text book examples, such as the electrostatic attraction/repulsion for a charge uniformly distributed on a long rod, or a large plate. Clearly the force does depend on the distribution of charge/mass.
OTOH, with regard to gravity, because gravity is such a weak force, most of the practical examples with rigid bodies in celestial mechanics do involve spherical bodies. One important exception is to treat the gravitational field of a spiral galaxy (it is not rigid, but its mass distribution can be treated as constant). This is not the same as the gravitational field of a central mass. I have shown how it can be treated in The effects of turbulence generated in Big Bang nucleosynthesis
A: The simple explanation is that any finite body (i.e. occupying a bounded region of space) looks like a point from sufficiently far away. This observation also tells you what is the range of validity of this "law". The distance between the bodies needs to be much larger than the linear size of each body. 
Using math and calculus it is possible to turn this intuition into precise and predictive equations. This approach goes under the name of multipole expansion
A: It is a fundamental "law" of nature that just happens to be true.
You can prove laws to a very high degree of accuracy, but from theory alone it is impossible to predict something fundamental by definition. Maths can't tell you the universe's axoims!
To give an example, Coulomb's law for the force between two charged masses is very similar to the gravitational law.
If this law is correct, it can be shown that the electric field inside any conductor is zero. Thus, instead of trying to test the law directly (separating to charges and measuring their forces very precisely), which is not very accurate, one can check the inverse square law for electrostatics to a very high precision by testing that the electric field really is zero inside a conductor by seeing if any charge flows onto a test rod.
Thus, although some laws may be hard to verify experimentally directly (such as by moving stars), we can check that the laws match up with our measurement of elliptical orbits and their time periods, amongst other things.
In fact, it was the orbit of Mercury that led scientists (Einstein) to show that Newton's Law of Gravitation was in fact incorrect as it did not match the experimental data.
So my answer to your question is that no, this law cannot be proven from some other law since it is fundamental, however, through experimental data, we can show that it is true and thus "good enough for the time being" until we find a better one! Can you really ever say that it is really true everywhere in the Universe without visiting every location in space and time? I'd say that's philosophy not Physics.  
