According to Newton's first principle there exists at least one
reference frame, called inertial, in which, by definition, if the
sum of all forces acting on a mass-point is zero, the point moves
uniformly along a straight line . Let us choose such a coordinate
frame. Any other coordinate frame that moves uniformly along a
straight line with respect to the given inertial reference frame
is also an inertial reference frame.
According to Newton's second postulate, in an inertial
reference frame the sum of all forces acting on a mass point
equals the acceleration of the point times its mass.
Newton's third principle states the if a mass point $A$ exerts
force $F$ on a mass point $B$ then, in return, $B$ exerts force
$-F$ on $A$ (recall that $F$ is a vector!).
We are given two mass points $A_1$ and $A_2$ with masses $m_1$ and
$m_2$ respectively. They generate a common gravitational field and
we want to study how the two points move with time. Let us fix an
inertial frame of reference $Oe_1e_2e_3$, which always exists
according to Newton's first principle. Here, $O$ is the origin of the coordinate system, while $e_1, e_2$ and $e_3$ are three unit vectors, perpendicular to each other. Let
$r_1=r_1(t)=x_1(t)e_1+y_1(t)e_2+z_1(t)e_3$ be the position vector of
point $A_1$ at time $t$. Analogously, let
$r_2=r_2(t)=x_2(t)e_1+y_2(t)e_2+z_2(t)e_3$ be the position vector of
$A_2$ at time $t$. According to Newton's law of gravity, mass
point $A_2$ attracts mass point $A_1$ with force
$$F_1(r_1,r_2)=
\left(\frac{m_1m_2G}{|r_2-r_1|^2}\right)\frac{r_2-r_1}{|r_2-r_1|}
=\frac{m_1m_2G}{|r_2-r_1|^3}(r_2-r_1)$$ where
$|r_2-r_1|=\sqrt{(x_2-x_1)^2+(y_2-y_1)^2+(z_2-z_1)^2}$ is the
distance between $A_1$ and $A_2$.
According to Newton's third principle, the force with which $A_1$
attracts $A_2$ is
$$F_2(r_1,r_2)=-F_1(r_1,r_2)=\frac{m_1m_2G}{|r_2-r_1|^3}(r_1-r_2)$$
By definition, for $j=1,2$ the velocity of the point $A_j$ is
$r^{\prime}_j = r'_j(t) = \frac{dr_j}{dt}$ and its acceleration is $r''_j
= r''_j(t) = \frac{d^2r_j}{dt^2}$. Consequently, Newton's second
law produces the following system of ordinary differential
equations
\begin{align*}
m_1 r''_1 &=F_1(r_1,r_2)\\
m_2 r''_2 &=F_2(r_1,r_2) = -F_1(r_1,r_2)
\end{align*}
or more explicitly
\begin{align*}
m_1 r''_1 &= \frac{m_1m_2G}{|r_2-r_1|^3}(r_2-r_1)\\
m_2 r''_2 &= \frac{m_1m_2G}{|r_2-r_1|^3}(r_1-r_2)
\end{align*}
Factor out the masses on both sides of the equations and obtain
\begin{align*}
r''_1 &= \frac{m_2G}{|r_2-r_1|^3}(r_2-r_1)\\
r''_2 &= \frac{m_1G}{|r_2-r_1|^3}(r_1-r_2)
\end{align*}
We end up with a system of two ordinary vector differential
equations or equivalently, a system of six scalar differential
equations. Here is where the magic happens and the system reduces
from six to three scalar equations describing the motion of one
point only.
Let us form the vector
$$R=\frac{m_1r_1+m_2r_2}{m_1+m_2}.$$ This vector $R=R(t)$
is the position of the common center of gravity of the two mass
points $A_1$ and $A_2$ at time $t$. After differentiating $R$
twice with respect to time $t$, we get
\begin{align*}
R^{''} &=\frac{1}{m_1+m_2}\Bigl(m_1r''_1+
m_2r''_2\Bigr)=\frac{1}{m_1+m_2}\Bigl(F_1+F_2\Bigr)\\
&=\frac{1}{m_1+m_2}\Big(F_1-F_1\Big)=0.
\end{align*}
which means that the velocity of the common center of mass of
points $A_1$ and $A_2$ is constant, i.e. $R^{'}(t) \equiv V_R$
which means that $R = R(t) = R_0 + tV_R$ moves uniformly along a
straight line. Then if we wish, according to Newton's first
principle, we can move the coordinate system to $R(t)$ and still
obtain an inertial reference frame. But let us keep working with
an arbitrary reference frame.
Now, let us set $r = r_2-r_1$. Then
\begin{align*}
r'' &= r_2''-r''_1 \\
&=-\frac{m_1G}{|r_2-r_1|^3}(r_2-r_1)-\frac{m_2G}{|r_2-r_1|^3}(r_2-r_1)\\
&=-\frac{(m_1+m_2)G}{|r_2-r_1|^3}(r_2-r_1)\\
&=-\left(\frac{(m_1+m_2)G}{|r|^3}\right)\,r
\end{align*}
Thus, we are left with a system of three scalar differential
equations instead of six.
$$\frac{d^2r}{dt^2}=-\left(\frac{(m_1+m_2)G}{|r|^3}\right)\,r$$
which can be written in terms of the reduced mass as
$$\left(\frac{m_1m_2}{m_1+m_2}\right)\frac{d^2r}{dt^2} = -\left(\frac{m_1m_2G}{|r|^3}\right)\,r$$
The solution $r=r(t)$ to this vector differential equation gives the solution to the
original problem
\begin{align*}
r_1(t) &=R(t) -\left(\frac{m_2}{m_1+m_2}\right)r(t)
=R_0+tV_R -\left(\frac{m_2}{m_1+m_2}\right)r(t)\\
r_2(t)&=R(t)+\left(\frac{m_1}{m_1+m_2}\right)r(t)
=R_0+tV_R+\left(\frac{m_1}{m_1+m_2}\right)r(t)
\end{align*}
Now, let us move our reference frame to the common center of
gravity of the two points $A_1$ and $A_2$. This means that $R_0 =
0$ and $V_R = 0$ and thus $R(t) \equiv 0$. Then the motion of the
two mass points becomes
\begin{align*}
r_1(t) &=-\left(\frac{m_2}{m_1+m_2}\right)r(t)\\
r_2(t)&=\left(\frac{m_1}{m_1+m_2}\right)r(t)
\end{align*}
If we re-express $r(t)$ in terms of $r_1(t)$ and $r_2(t)$ we obtain
$$r(t) = - \left(\frac{m_1 + m_2}{m_2}\right)r_1(t) = \left(\frac{m_1+m_2}{m_1}\right)r_2(t)$$
and thus, the system of differential equations in the reference
frame with respect to the common center of gravity of the two mass
points becomes (by plugging each of these two expression for $r$
into the system for $r$)
\begin{align*}
r''_1&= - \left(\frac{(m_2)^3G}{(m_1+m_2)^2|r_1|^3}\right)r_1\\
r''_2&=-\left(\frac{(m_1)^3G}{(m_1+m_2)^2|r_2|^3}\right)r_2
\end{align*}
which as you can see is already a decoupled system of two vector equations.
If the mass $m_2$ is much bigger than $m_1$, so that $m_1$ can be ignored, then $(m_1 + m_2) \sim m_2$, and thus
\begin{align*}
m_1 r''_1&= - \left(\frac{m_1m_2G}{|r_1|^3}\right)r_1\\
\end{align*}
is the motion of a mass point in a gravitational field created by a much bigger mass. This is what Galileo observed in his experiments -- the motion of a body in Earth's gravity does not depend on its mass. Well, technically, Galileo was wrong, as we can see from our equations, but for practical purposes he was right. He didn't have precise enough measurements to see the difference between $(m_1+m_2)$ and $m_1$.
How do we derive
\begin{align*}
r_1(t) &=R(t) - \left(\frac{m_2}{m_1+m_2}\right)r(t)
=R_0+tV_R+\left(\frac{m_2}{m_1+m_2}\right)r(t)\\
r_2(t)&=R(t) + \left(\frac{m_1}{m_1+m_2}\right)r(t)
=R_0+tV_R-\left(\frac{m_1}{m_1+m_2}\right)r(t)?
\end{align*}
Think of vector $R$ as a three dimensional arrow with origin $O$ -- which is the origin of the coordinate system $Oe_1e_2e_3$ -- and end the center of mass of segment $A_1A_2$. This is called a position vector of the center of mass. The vector $r_j$ is the analogous position vector of point $A_j$ -- starts from $O$ and ends at $A_j$, for $j=1,2$. Now the relative position vector $r = r_2-r_1$ is the vector with origin $A_1$ and end $A_2$. The line determined by the mass points $A_1$ and $A_2$ is $A_1A_2$, which is the line passing through the center of gravity $R$ and with vector $r$ lying on it, and its parametric equation is $p(\lambda) = R + \lambda r$. Therefore $r_j = R + \lambda_j r$ for some number $\lambda_j$ where $j=1,2$. We have to determine $\lambda_1$ and $\lambda_2$. Look at the systems
\begin{align}
r_1 &= R + \lambda_1 r\\
r_2 &= R + \lambda_2 r
\end{align}
and
\begin{align}
R &= R=\frac{m_1r_1+m_2r_2}{m_1+m_2}\\
r &= r_2 - r_1
\end{align}
Plug the first system into the second system and you find equations for the $\lambda$'s. when you solve them you get the expression above.
Let us come back to the system that features the reduced mass
$$\left(\frac{m_1m_2}{m_1+m_2}\right)\frac{d^2r}{dt^2} = -\left(\frac{m_1m_2G}{|r|^3}\right)\,r$$
which I am going to rewrite as
$$\mu \frac{d^2r}{dt^2} = -\kappa(r)\,r$$
where $\kappa(r) = \frac{m_1m_2G}{|r|^3}$ is a scalar function and $\mu = \frac{m_1m_2}{m_1+m_2}$. Let $$K(r') = \frac{1}{2} \mu |r'|^2 \,\,\text{ and }\,\, U(r) = -\frac{m_1m_2G}{|r|}$$ be the kinetic and the potential energies respectively. Now, you can form the total energy
$$H(r,r') = K(r') + U(r)$$ as well as the Lagrangian
$$\mathcal{L}(r,r') = K(r') - U(r).$$
The momentum of a mass point is usually linked to the Lagrangian by $p = \frac{\partial \mathcal{L}}{\partial\, r'} = \mu \, r'$. So indeed, in our case $p = \mu r'$. A direct computation shows that $H(r,r')$ (technically it should be $H(r,p)$) is a constant of motion, and it is in fact the Hamiltonian of our reduced system. Furthermore, one can form the angular momentum $L = r \times p = r \times (\mu r')$ so if we differentiate it with respect to $t$ we get
$$L' = r' \times (\mu r') + r \times (\mu r'') = r \times (\mu r'') = r \times (-\kappa(r) r) = -\kappa(r) (r \times r) = 0.$$ Hence, the angular momentum is also conserved. That basically tells you that the motion is planar. And you have conservation of energy. And if you look at the Runge-Lenz vector
$$A = P \times L - \mu \, m_1m_2\, G \, \frac{r}{|r|}$$ this is also conserved, which is basically the fact that the major axis of our motion is also constant with time (for elliptical orbit, this is the major axis of the ellipse. It stays unchanged with time). So everything goes down as you have seen it. You can fix your plane of motion, determined by the constant angular momentum $L$, use the Runge-Lenz vector as an axis for your polar coordinates in that plane, and finally use the conservation of energy to derive your simple equation of motion in polar coordinates, where the mass will be the reduced mass. So, your angle $\theta(t)$ is the angle between the motion vector $r = r_2-r_1$ and the Runge-Lenz vector, all moving in the plane of constant angular momentum. Your radial coordinate $\rho(t) = |r(t)| = |r_2(t)-r_1(t)|$ is the distance between the two mass points. These coordinates are not a direct switch from Cartesian inertial frame of reference to some polar reduction of it. They are very closely related to the inertial frame, but the center of the polar coordinates is not exactly the center of mass. You think of them as slightly abstract coordinates, although, as you can see, they have a clear geometric meaning. If you want a direct interpretation simply take one of the equations
$$r^{''}_1 = - \left(\frac{(m_2)^3G}{(m_1+m_2)^2|r_1|^3}\right)r_1$$ and apply absolutely the same arguments as before (total energy, angular momentum, Runge-Lenz vector), and you get your polar coordinate to be exactly a direct switch from Cartesian inertial frame to a polar restriction on the constant angular momentum plane. Here, the center of the polar system is exactly the center of mass (coinciding with the focus of your orbit, which is a conic section), the Runge-Lenz vector is the zero axis of your polar coordinates. But the reduced mass is mixed up in the constants somewhere. Whichever you prefer.
Why does not the precession rate of the angular velocity vector give
exactly the frequency of wobbling? In other words, how come the
precession of $\vec\omega$ is different to the line of nodes rotation?
Short answer: Because the precession rate $\Omega$ corresponds to the precession of $\vec\omega$ in the body frame, not in the Earth frame. The variables in the Euler equations, $I_i$, $\omega_i$ and $\tau_i$, are in the body frame.
The angular velocity of the disk has two contributions, the component $\vec\omega'$ due to the spin and the component $\vec\omega''$ due to a rotation of the inclined disc about a vertical axis. The latter corresponds to the wobbling as viewed by someone in the Earth frame.
As we can see from the figure, the resultant angular velocity $\vec\omega$ is always off the symmetry axis of the disk, which means it has a non-vanishing projection in the plane of the disk (dashed line). At same time, the axes fixed in the plane of the disk, $\rho_1$ and $\rho_2$ rotates with spin $s$ (viewed by someone in the Earth frame). Hence, someone in the disk frame will see the axes fixed and the projection of $\omega$ in this plane rotating with rate $s$. That is why the precession rate equals the spin.
On the other hand, the wobbling effect is due to the non-vanishing $\vec\omega''$. A $2\pi$ rotation about the vertical line corresponds to a complete oscillation (wobbling) therefore the frequency of the wobbling is actually equal to the magnitude of $\vec\omega''$.
Best Answer
ja72's answer is probably right, but I was confused by his notation, so I will give my own answer. Suppose we have two object rotating with angular velocity $\vec{\omega}_1$ and $\vec{\omega}_2$. Then the velocity of a point $\vec{r}$ of object $1$ in the lab frame is $\vec{v}_{1,lab}=\vec{\omega}_1 \times \vec{r}$. Similarly, the velocity of a point $\vec{r}$ of object $2$ in the lab frame is $\vec{v}_{2,lab}=\vec{\omega}_2 \times \vec{r}$.
Now to someone in the second object, a point $\vec{r}$ that is stationary in the lab frame will have an apparent velocity $-\vec{v}_{2,lab}=-\vec{\omega}_2 \times \vec{r}$. From this, you can see that a point $\vec{r}$ in the first object will appear to have a velocity $\vec{v}_{1,lab}-\vec{v}_{2,lab} = \vec{\omega}_1 \times \vec{r}-\vec{\omega}_2 \times \vec{r} = (\vec{\omega}_1-\vec{\omega}_2)\times \vec{r}$. Thus the first object appears to have angular velocity $\vec{\omega}_1-\vec{\omega}_2$ to an observer in the second object.
Now since $\vec{\omega}_1$ is stationary in the lab frame, it's time derivative is $- \vec{\omega}_2 \times\vec{\omega}_1 =\vec{\omega}_1\times \vec{\omega}_2$ in the frame of the second object. Thus object one appears to have an angular acceleration of $\vec{\omega}_1\times \vec{\omega}_2$ in the frame of the second object.