Q:

Euler-Lagrange equations for free rigid body

How do I derive Euler's equations of motion for a free rigid body using a Lagrangian formulation? 

The required equations are, in vector form, $J\dot{\omega} = -\omega \times J\omega$ where $J$ is moment of inertia of body and $\omega$ is angular velocity
My attempt:
$L(\omega) = \frac{1}{2}\omega^TJ\omega$ is the Lagrangian.
Using Euler Lagrange equation, $\frac{d}{dt}\frac{\partial L}{\partial \dot{q}} - \frac{\partial L}{\partial q} = 0 $,
$\dot{q} = \omega$
$\frac{\partial L}{\partial \omega} = J\omega$
$\frac{d}{dt}\frac{\partial L}{\partial \dot{q}} = \frac{d}{dt}\frac{\partial L}{\partial \omega} = \frac{d}{dt}(J\omega)$
$\frac{\partial L}{\partial q} = 0 $
So finally, $\frac{d}{dt}(J\omega) = 0 \implies \omega = 0 $ , but this is not the correct equation. 
I require $\frac{d}{dt}\frac{\partial L}{\partial \omega} = \frac{\partial L}{\partial \omega} \times \omega$ as the final equation which is equivalent to the required equation since $\frac{\partial L}{\partial \omega} = J\omega$. 

A:

In general $J$ is a function of time when expressed in coordinates fixed to the inertial frame (as it is in your derivation- and same with $\omega$). This is easy to see: $J$ describes the distribution of the body's mass about the inertial coordinates, so as the body rotates relative to the inertial coordinates, this distribution will change. However, since the body is rigid, $J$ expressed in coordinates fixed to the body itself should be constant. Let $R$ be the rotation matrix taking body coordinates to inertial coordinates, and let $J_B$ be $J$ expressed in body coordinates. Then,
$$
\begin{align}
\frac{d}{dt}(J\omega) &= J\dot{\omega} + \dot{J}\omega\\
&= J\dot{\omega} + \dot{(R J_B R^\intercal)}\omega\\
&= J\dot{\omega} + \dot{R} J_B R^\intercal\omega + R \dot{J_B} R^\intercal\omega + R J_B \dot{R^\intercal}\omega\\
&= J\dot{\omega} + \Omega R J_B R^\intercal\omega - R J_B R^\intercal \Omega\omega\\
&= J\dot{\omega} + \omega \times J\omega
\end{align}
$$
since $\dot{R} = \Omega R$ where $\Omega$ is the skew-symmetric matrix of $\omega$ (i.e. its cross product operator), $\Omega \omega = 0$, and $\dot{J_B}=0$. Since your Lagrangian showed that $\frac{d}{dt}(J\omega) = 0$, then we have, as you expected,
$$
J\dot{\omega} = -\omega \times J\omega
$$
In general, we let $\tau$ be the net external torque on the body and say,
$$
\frac{d}{dt}(J\omega) = \tau
$$
It is interesting that this equation doesn't change form even if we change our basis for all vectors to body coordinates by left multiplying $R^\intercal$ and defining $\omega = R \omega_B$ and $\tau = R \tau_B$,
$$
\begin{align*}
R^\intercal\big{(}J\dot{\omega} + \omega \times J\omega\big{)} &= R^\intercal \tau\\
R^\intercal (R J_B R^\intercal) (R \dot{\omega}_B) + R^\intercal (R J_B R^\intercal) (\dot{R} \omega_B) + R^\intercal \big{(}(R \omega_B) \times (R J_B R^\intercal) (R \omega_B)\big{)} &= \tau_B\\
J_B \dot{\omega}_B + J_B R^\intercal \Omega \omega + \omega_B \times J_B \omega_B &= \tau_B\\
J_B \dot{\omega}_B + \omega_B \times J_B \omega_B &= \tau_B
\end{align*}
$$
Note that $R$ can only "factor out" of a cross-product because it is orthogonal.
Lastly, to be fully convinced of the $\dot{R}$ identity, just take a time-derivative of the orthogonality identity $R^\intercal R = I$. In body coordinates we have $\dot{R} = R \Omega_B$.

