Great question; I remember being so confused by this when I first took analytic mechanics.
The components of the angular velocity "in the body frame" aren't zero because when one writes these components, one isn't referring to measurements of the motions of the particles in the body frame (because, of course, the particles are stationary in this frame). Instead, one is referring to angular velocity as measured in an inertial frame but whose components have simply been written with respect to a time-varying basis that is rotating with the body.
In practice, we make measurements of the positions $\mathbf x_i(t)= (x_i(t),y_i(t),z_i(t))$ of the particles in an inertial frame. Then, we note that for a rigid body (let's consider pure rotation for simplicity), the position of each particle $i$ satisfies
\begin{align}
\mathbf x_i(t) = R(t) \mathbf x_i(0)
\end{align}
for some time-dependent rotation $R(t)$. Then we compute $\boldsymbol\omega(t) = (\omega^x(t),\omega^y(t),\omega^z(t))$ in the standard way in terms of $R(t)$. To see how this is done in detail, see, for example
https://physics.stackexchange.com/a/74014/19976
Once we have $\boldsymbol\omega$, we can write its components with respect to any basis we like. If we write it in the standard ordered basis $\{\mathbf e_i\}$, then we'll just get $\omega_x(t)$ as its components. If we write it in some basis $\{\mathbf e_{i,B}(t)\}$ that is rotating with the body (like one that points along the principal axes of the body) then we get different components $\omega^i_B(t)$, and these are the body components.
Main Point Reiterated. Angular velocity is being measured with respect to an inertial frame, but its components can be taken with respect to any basis we wish such as one rotating with the body.
Do I need to use the angular velocity vector in the rotating or inertial reference frame for this?
Yes. You can do it either way.
I start with the expression that relates the time derivative of a vector quantity $\boldsymbol u$ in the inertial and rotating frames:
$$\left(\frac {d\boldsymbol u}{dt}\right)_I
= \left(\frac {d\boldsymbol u}{dt}\right)_R
+ \boldsymbol \omega \times \boldsymbol u$$
Expressing this in coordinate spaces,
$$\dot{\boldsymbol u}_I
= \mathrm T_{R\to I}
\left(\dot{\boldsymbol u}_R
+ \boldsymbol \omega_R \times \boldsymbol u_R
\right)
= \mathrm T_{R\to I}\,\dot{\boldsymbol u}_R
+ \mathrm T_{R\to I} \mathrm{Sk}[\boldsymbol \omega_R]\,\boldsymbol u_R
$$
where $\mathrm T_{R\to I}$ is the matrix that transforms a vector from its representation in the rotating frame to its corresponding representation in the inertial frame and $\mathrm{Sk}[\boldsymbol x]$ is the skew symmetric such that $\mathrm{Sk}[\boldsymbol x]\,\boldsymbol a = \boldsymbol x \times \boldsymbol a$ (in other words, the cross product matrix).
Differentiating $\boldsymbol u_I = \mathrm T_{R\to I}\,\boldsymbol u_R$ with respect to time yields
$$
\dot{\boldsymbol u}_I
= \mathrm T_{R\to I}\,\dot{\boldsymbol u}_R
+ \dot{\mathrm T} _{R\to I}\,\boldsymbol u_R
$$
From which $\dot{\mathrm T} _{R\to I} = \mathrm T_{R\to I}\,\mathrm{Sk}[\boldsymbol \omega_R]$. You can similarly derive that $\dot{\mathrm T} _{I\to R} = -\mathrm T_{I \to R}\,\mathrm{Sk}[\boldsymbol \omega_I]$. Transposing yields alternate expressions. Summarizing,
$$\begin{aligned}
\dot{\mathrm T} _{R\to I}
&= \mathrm T_{R\to I}\,\mathrm{Sk}[\boldsymbol \omega_R] \\
&= \mathrm{Sk}[\boldsymbol \omega_I]\,\mathrm T_{R\to I} \\
\dot{\mathrm T} _{I \to R}
&= - \mathrm T_{I \to R}\,\mathrm{Sk}[\boldsymbol \omega_I] \\
&= - \mathrm{Sk}[\boldsymbol \omega_R]\,\mathrm T_{I \to R}
\end{aligned}$$
When numerically integrating this, together with Euler's equation of rotation, is there a way to ensure that the determinant of R remains equal to one (otherwise x⃗ (t) will also be scaled)?
Yes, using Lie group methods. Explaining this is a bit long. 140 pages long. I'll instead refer you to Iserles, Arieh, et al. "Lie-group methods." Acta Numerica 2000 9 (2000): 215-365.
Best Answer
Yes. It's one of many consequences of Euler's rotation theorem. This neat little trick of finding the Euler rotation axis works in three dimensional space, and in three dimensional space only. The Lie group SO(3) is a very special space.
A more generic way of looking at rotation is that rotation in any Euclidean space comprises a combinations of rotations parallel to a two dimensional plane rather than about an axis. There's only one plane in two dimensional space, so rotation in 2D space requires but one parameter. There are three pairs of planes of rotation in three dimensional space (the YZ, ZX, and XY planes), so rotation in 3D space requires three parameters. You can think of those planar rotations as being about an axis (with an arbitrary decision as to what qualifies as positive or negative rotation). There are six planes of rotation in four dimensional space, so rotations in 4D space requires six parameters. This can result in something very weird, something you don't see in 2D or 3D space, a Clifford rotation.
Source: http://eusebeia.dyndns.org/4d/vis/10-rot-1