Lagrangian equations of motion for ball rolling on turntable The equations governing the motion of a ball of mass $m$, radius $R$ rolling on a table rotating at constant angular velocity $ \Omega $ which are derived using Newton's laws are: (I present these for comparison) 
\begin{align*}
(I+mR^2) \dot{\omega}_x &= (I+2mR^2)\Omega \omega_y-mR \Omega^2y \\
(I+mR^2) \dot{\omega}_y &= -(I+2mR^2)\Omega \omega_x + mR \Omega^2 x \\
\dot{x} &= R \omega_y \\
\dot{y} &= -R \omega_x
\end{align*}
Where $x,y,\omega_x,\omega_y$ are absolute values measured in the rotating frame ($x,y$ being positions and $\omega_x,\omega_y$ angular velocities of the ball). To express the position, velocity, etc in the inertial $XYZ$ frame we can perform a change of variables:
\begin{align*}
X=x \cos \theta - y \sin \theta \\
Y=x \sin \theta + y \cos \theta
\end{align*}
The equations written above are correct as far as I know. Now I've tried to derive these equations using the Lagrangian approach, but my equations differ slightly from the above. I'll share my work here:
We start with the standard formulation of the Lagrangian equations of motion:
$$\frac{d}{dt} \left( \frac{\partial L}{\partial \dot{q_i}} \right) - \frac{\partial L}{\partial q_i}=Q_i$$
For this system, there are no non-conservative forces doing work at any time so $Q_i=0$ (Assuming no-slip here). Now the kinetic energy of the system is:
$$T=\frac{1}{2}m \| v \| ^2 + \frac{1}{2} I \| \omega \|^2  $$
Where $v$,$\omega$ are the absolute linear and angular velocities of the ball, $I$ is the moment of inertia of the center of mass. I'll proceed in the rotating frame basis $xyz$, where the position $r$ of the ball is given by
$$\vec{r}=x \hat{i} + y \hat{j}$$
Using velocity kinematics in rotating frames, the absolute velocity of the ball is given by:
$$\vec{v}=\vec{\Omega} \times \vec{r} + \vec{v}_{xyz}$$
Where $\vec{v}_{xyz}$ is the velocity of the ball in the rotating frame:
$$\vec{v}_{xyz} = \dot{x} \hat{i} + \dot{y} \hat{j}$$
So the absolute velocity $\vec{v}$ after expanding gives:
$$\vec{v} = (\dot{x}-\Omega y) \hat{i} + (\dot{y} +\Omega x) \hat{j} $$
It follows from taking the magnitude of the absolute velocity:
$$\| v \|^2 = (\dot{x}-\Omega y)^2 + (\dot{y}+\Omega x)^2$$
The first unknown term in the Lagrangian equations of motion. The absolute angular velocity $\vec{\omega}$ is more straight-forward:
$$\vec{\omega} = \vec{\Omega} + \vec{\omega}_{xyz}$$
Where $\vec{\omega}_{xyz}$ is the angular velocity of the ball in the rotating frame, where one can show it DOES NOT have a component in the $\hat{k}$ direction if we impose no-slip kinematics (which we are), so:
$$\vec{\omega}_{xyz} = \omega_x \hat{i} + \omega_y \hat{j}$$
And since 
$$\vec{\Omega} = \Omega \hat{k}$$
The absolute angular velocity is:
$$\vec{\omega} = \omega_x \hat{i} + \omega_y \hat{j} + \Omega \hat{k} $$
Then it follows that:
$$\| \omega \|^2 = \omega_x^2 + \omega_y^2 + \Omega^2$$
The potential energy of the system is constant and doesn't affect the equations of motion so, the Lagrangian becomes:
$$L=\frac{1}{2} m \left[ (\dot{x}-\Omega y)^2 + (\dot{y}+\Omega x)^2 \right] + \frac{1}{2} I \left[ \omega_x^2 + \omega_y^2 + \Omega^2 \right] $$
Now the constraint equations are the no-slip conditions:
$$\vec{v}_{xyz} = \vec{\omega}_{xyz} \times \vec{R}$$
Where $\vec{R} = R \hat{k} $ of course, then we have two conditions:
\begin{align*}
\dot{x} &= R \omega_y \\
\dot{y} &= - R \omega_x
\end{align*}
Which are non-holomonic constraints, but given their simple nature I opted out of doing Lagrangian multipliers and simply substituted them into the equations of motion (I reworked it with Lagrangian multipliers after and got the same thing). Now after substituting the constraints, the Lagrangian becomes:
$$L=\frac{1}{2} \left( m + \frac{I}{R^2} \right) (\dot{x}^2 + \dot{y}^2) - m \Omega \dot{x} y + m \Omega \dot{y} x + \frac{1}{2} m \Omega^2 (x^2 + y^2) + \frac{1}{2} I \Omega^2 $$
From here, applying the Lagrangian equation above with $Q_i=0$ I get the following equations of motion:
$$\left( m+\frac{I}{R^2} \right) \ddot{x} - 2m \Omega \dot{y} - m \Omega^2 x = 0 $$
$$\left( m+\frac{I}{R^2} \right) \ddot{y} - 2m \Omega \dot{x} + m \Omega^2 y = 0$$
And using the no-slip conditions again we can rewrite:
\begin{align*}
\left( I + mR^2 \right) \dot{\omega}_x &= 2mR^2 \Omega \omega_y - mR \Omega^2 y \\
\left( I + mR^2 \right) \dot{\omega}_y &= -2mR^2 \Omega \omega_x + mR \Omega^2 x 
\end{align*}
Now if you compare these last two equations with the ones I wrote in the beginning, the only difference is in the first term on the right-hand side. Look at these two for instance:
\begin{align*}
(I+mR^2) \dot{\omega}_x &= (I+2mR^2)\Omega \omega_y-mR \Omega^2y \\
\left( I + mR^2 \right) \dot{\omega}_x &= 2mR^2 \Omega \omega_y - mR \Omega^2 y
\end{align*}
The ONLY difference is the missing $I$ term! I'm missing the moment of inertia for some reason, why is that? What's wrong about my Lagrangian approach?


Newton Euler Equations
\begin{align*}
& \underbrace{\left[ \begin {array}{cccc} m&0&0&0\\  0&m&0&0
\\  0&0&{\it I_s}&0\\  0&0&0&{\it I_s}
\end {array} \right]}_{\boldsymbol A}
\begin{bmatrix}
  \ddot{x} \\
  \ddot{y} \\
  \dot{\omega_x} \\
  \dot{\omega_y} \\
\end{bmatrix}=
\underbrace{\begin{bmatrix}
              -m \left( -{\Omega}^{2}x-2\,\Omega\,{\dot y} \right)  \\
              -m \left( -{\Omega}^{2}y+2\,\Omega\,{\dot x} \right)  \\
              0 \\
              0 \\
            \end{bmatrix}
}_{\boldsymbol b}+\boldsymbol C_B^T\,\boldsymbol\lambda_n\\\\
\end{align*}
where $~I_s=m\,\frac 25\,\rho^2~$  the sphere moment of inertia and $\rho~$  the sphere radius, $~\boldsymbol C_B^T~$ distribution matrix of the generalized constraint forces $~ \boldsymbol\lambda_n$
Kinematic
\begin{align*}
 &\begin{bmatrix}
  \dot{x} \\
  \dot{y} \\
  {\omega_x} \\
  {\omega_y} \\
\end{bmatrix}=\underbrace{\left[ \begin {array}{cc} 0&-\rho\\ \rho&0
\\  1&0\\  0&1\end {array} \right]}_{\boldsymbol J}
\begin{bmatrix}
  \omega_x\\
  \omega_y \\
\end{bmatrix}
\end{align*}
The equations of motion
\begin{align*}
&\boldsymbol J^T\,\boldsymbol A\,\boldsymbol J\,\begin{bmatrix}
  \dot\omega_x\\
  \dot\omega_y \\
\end{bmatrix}=\boldsymbol J^T\,\boldsymbol b\\
&\Rightarrow\\\\
&\begin{bmatrix}
     \dot{\omega}_x \\
     \dot{\omega}_y \\
   \end{bmatrix}=\left[ \begin {array}{c} {\frac {\rho\,m\Omega\, \left( -\Omega\,y+2
\,\omega_{{y}}\rho \right) }{{\rho}^{2}m+I{{s}}}}
\\ -{\frac {\rho\,m\Omega\, \left( -\Omega\,x+2\,
\omega_{{x}}\rho \right) }{{\rho}^{2}m+I_{{s}}}}\end {array} \right]
\tag{1}\\
&\text{and}\\
&\dot{x}=-\rho\,\omega_y  \tag{2}\\
&\dot{y}=\rho\,\omega_x \tag{3} \\
\end{align*}
Initial system Sphere Position
\begin{align*}
 & \boldsymbol R_s=\left[ \begin {array}{ccc} \cos \left( \Omega\,t \right) &-\sin
 \left( \Omega\,t \right) &0\\ \sin \left( \Omega\,t
 \right) &\cos \left( \Omega\,t \right) &0\\ 0&0&1
\end {array} \right]
\begin{bmatrix}
     x \\
     y \\
     \rho \\
   \end{bmatrix}
\end{align*}
the  rotation matrix $\boldsymbol S~$ between local coordinate system and inertial system is:
\begin{align*}
  &\boldsymbol S=\left[ \begin {array}{ccc} \cos \left( \varphi _{{z}} \right) &-\sin
 \left( \varphi _{{z}} \right) &0\\  \sin \left(
\varphi _{{z}} \right) &\cos \left( \varphi _{{z}} \right) &0
\\  0&0&1\end {array} \right]
\, \left[ \begin {array}{ccc} \cos \left( \varphi _{{y}} \right) &-\sin
 \left( \varphi _{{y}} \right) &0\\  \sin \left(
\varphi _{{y}} \right) &\cos \left( \varphi _{{y}} \right) &0
\\  0&0&1\end {array} \right]
\,\left[ \begin {array}{ccc} \cos \left( \varphi _{{x}} \right) &-\sin
 \left( \varphi _{{x}} \right) &0\\  \sin \left(
\varphi _{{x}} \right) &\cos \left( \varphi _{{x}} \right) &0
\\  0&0&1\end {array} \right]\\\
&\text{with}\\
&\boldsymbol{\dot{S}}=\left[ \begin {array}{ccc} 0&-\omega_{{z}}&\omega_{{y}}
\\  \omega_{{z}}&0&-\omega_{{x}}\\
-\omega_{{y}}&\omega_{{x}}&0\end {array} \right]
\boldsymbol S\,\\\\
&\Rightarrow\\
&\begin{bmatrix}
   \omega_x \\
   \omega_y \\
   \omega_z \\
 \end{bmatrix}=\left[ \begin {array}{ccc} \cos \left( \varphi _{{z}} \right) \cos
 \left( \varphi _{{y}} \right) &-\sin \left( \varphi _{{z}} \right) &0
\\  \sin \left( \varphi _{{z}} \right) \cos \left(
\varphi _{{y}} \right) &\cos \left( \varphi _{{z}} \right) &0
\\  -\sin \left( \varphi _{{y}} \right) &0&1
\end {array} \right]
\,\begin{bmatrix}
    \dot\varphi_x \\
    \dot\varphi_y \\
    \dot\varphi_z \\
  \end{bmatrix}
 \end{align*}
from here you obtain
\begin{align*}
  \begin{bmatrix}
    \dot{\varphi_x}\\
    \dot{\varphi_y} \\
    \dot{\varphi_z} \\
  \end{bmatrix}=
  \left[ \begin {array}{ccc} {\frac {\cos \left( \varphi _{{z}}
 \right) }{\cos \left( \varphi _{{y}} \right) }}&{\frac {\sin \left(
\varphi _{{z}} \right) }{\cos \left( \varphi _{{y}} \right) }}&0
\\  -\sin \left( \varphi _{{z}} \right) &\cos \left(
\varphi _{{z}} \right) &0\\  {\frac {\cos \left(
\varphi _{{z}} \right) \sin \left( \varphi _{{y}} \right) }{\cos
 \left( \varphi _{{y}} \right) }}&{\frac {\sin \left( \varphi _{{z}}
 \right) \sin \left( \varphi _{{y}} \right) }{\cos \left( \varphi _{{y
}} \right) }}&1\end {array} \right]
\begin{bmatrix}
  \omega_x \\
  \omega_y \\
  \Omega \\
\end{bmatrix}
 \tag {4}
 \end{align*}
Altogether you obtain  7 first order differential equations
Euler Lagrange with non holonomic constraint equations
\begin{align*}
   &\frac{d}{dt}\left(\frac{\partial \mathcal{L}}{\partial \dot{\boldsymbol{w}}}\right)^T-\left(
  \frac{\partial \mathcal{L}}{\partial \boldsymbol{w}}\right)^T=\left[\frac{\partial \boldsymbol{R}}{\partial \boldsymbol{w}}\right]^T\,\boldsymbol{F}_s+
   \left(\frac{\partial \boldsymbol{g}_n}{\partial \dot{\boldsymbol{w}}}\right)^T\boldsymbol{\lambda}_n\tag A
 \end{align*}
where
\begin{align*}
& \boldsymbol{R}=\begin{bmatrix}
                  x \\
                  y \\
                \end{bmatrix}\\
 & \dot{\boldsymbol{w}}=\begin{bmatrix}
                \dot{x} \\
                \dot{y} \\
                \omega_x \\
                \omega_y \\
              \end{bmatrix}~,
              {\boldsymbol{w}}=\begin{bmatrix}
                {x} \\
                {y} \\
                \dot{\varphi}_x \\
                \dot{\varphi}_y \\
              \end{bmatrix}\\\\
 &\mathcal{L}=\frac{m}{2}\left(\dot{x}^2+\dot{y}^2\right)+
 \frac{I_S}{2}\left[\omega_x~,\omega_y~,\Omega\right]^T\,
 \left[\omega_x~,\omega_y~,\Omega\right]  \\\\
 &\boldsymbol F_s=\begin{bmatrix}
                   -m \left( -{\Omega}^{2}x-2\,\Omega\,{\dot y} \right)  \\
              -m \left( -{\Omega}^{2}y+2\,\Omega\,{\dot x} \right)  \\
                  \end{bmatrix}
  \\\\
&\text{and the non holonomic constraint equations }\\
&\boldsymbol g_n=  \left[ \begin {array}{c} {\dot{y}}-\rho\,\omega_{{x}}
\\  {\dot{x}}+\rho\,\omega_{{y}}\end {array} \right]
 \end{align*}
from equation (A) you obtain:
\begin{align*}
  &\underbrace{\left[ \begin {array}{cccccc} m&0&0&0&0&-1\\ 0&m&0&0
&-1&0\\  0&0&I_s&0&\rho&0\\  0&0
&0&I_s&0&-\rho\\  0&-1&\rho&0&0&0
\\  -1&0&0&-\rho&0&0\end {array} \right]}_{\boldsymbol A_L}
\underbrace{\begin{bmatrix}
                \ddot{x} \\
                \ddot{y} \\
                \dot\omega_x \\
                \dot\omega_y \\
                \lambda_n
              \end{bmatrix}}_{\boldsymbol{\ddot{w}}}
=\underbrace{\begin{bmatrix}
               -m \left( -{\Omega}^{2}x-2\,\Omega\,{\dot y} \right)  \\
              -m \left( -{\Omega}^{2}y+2\,\Omega\,{\dot x} \right)  \\
               0 \\
               0 \\
               0 \\
             \end{bmatrix}
 }_{\boldsymbol b_L}\tag{B}  \\\
&\Rightarrow
 \end{align*}
substitute $~\dot{y}=\rho\,\omega_x~,\dot{x}=-\rho\,\omega_y~$ into equation (B) and solve for $~\ddot{\boldsymbol{w}}~$ you obtain
\begin{align*}
   &\begin{bmatrix}
     \dot{\omega}_x \\
     \dot{\omega}_y \\
   \end{bmatrix}=\left[ \begin {array}{c} {\frac {\rho\,m\Omega\, \left( -\Omega\,y+2
\,\omega_{{y}}\rho \right) }{{\rho}^{2}m+I{{s}}}}
\\ -{\frac {\rho\,m\Omega\, \left( -\Omega\,x+2\,
\omega_{{x}}\rho \right) }{{\rho}^{2}m+I_{{s}}}}\end {array} \right]
 \end{align*}
those  results are equal to  the results of equation (1) .
Edit
\begin{align*}
&\textbf{how to obtain the Jacobi-Matrix $~\boldsymbol J~$ (kinematic equations) }\\\\
 &\text{the non holonomic constraint equations are}\\
 &\boldsymbol g_n=\begin{bmatrix}
           \dot{x}-\rho\omega_y \\
           \dot{y}+\rho\omega_x \\
         \end{bmatrix}=\boldsymbol 0\\\\
&\text{obtain the time derivative }\\
&\boldsymbol{\dot{q}}_n=
 \underbrace{\left[ \begin {array}{cccc} 1&0&0&-\rho\\  0&1&\rho&0
\end {array} \right]}_{\boldsymbol C_n}
\underbrace{\begin{bmatrix}
  \dot{x} \\
  \dot{y} \\
  \omega_x \\
  \omega_y \\
\end{bmatrix}}_{\boldsymbol{\dot{w}}}=\boldsymbol C_n\,\boldsymbol{\dot{w}}=\boldsymbol 0\tag{c}
\end{align*}
you have two constraint equations and four velocities $\boldsymbol{\dot{w}}~$ thus the generalized velocities  are two (4-2) if you chose the generalized velocities $~\omega_x~,\omega_y$ you can obtain from equation (c)
\begin{align*}
&\begin{bmatrix}
  \dot{x} \\
  \dot{y} \\
  \omega_x \\
  \omega_y \\
\end{bmatrix}= \underbrace{\left[ \begin {array}{cc} 0&\rho\\ -\rho&0
\\  1&0\\  0&1\end {array} \right]}_{\boldsymbol J}\begin{bmatrix}
         \omega_x \\
         \omega_y \\
       \end{bmatrix}\\\
&\text{according to d'Alemberts principal }\\
&\boldsymbol J^T\,\boldsymbol C^T_n=
\left[ \begin {array}{cccc} 0&-\rho&1&0\\ \rho&0&0&
1\end {array} \right]
\,\left[ \begin {array}{cc} 1&0\\  0&1
\\  0&\rho\\  -\rho&0\end {array}
 \right]= \left[ \begin {array}{cc} 0&0\\ 0&0\end {array}
 \right]     ~  \surd
\end{align*}
