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$.

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 a rotation matrix taking body coordinates to inertial coordinates, and let $J_B$ be $J$ expressed body coordinates. Then, \begin{align} \frac{d}{dt}(J\omega) &= J\dot{\omega} + \dot{J}\omega\\ &= J\dot{\omega} + \dot{(R J_B R^T)}\omega\\ &= J\dot{\omega} + \dot{R} J_B R^T\omega + R \dot{J_B} R^T\omega + R J_B \dot{R^T}\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), 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$$