Non-Holonomic constraint in rigid body dynamics I have solved many problems on Holonomic constraint using Lagrange multiplier method but I don't know how to tackle problems on non-Holonomic constraint.
Can anyone help me with the following problem which was used my professor to illustrate non-Holonomic constraint in the lecture.

"A uniform rod of length $L$ and mass $M$, moving in a plane say $xy$. Now at one end it has a knife edge constraint which prevents velocity components perpendicular to the rod at that point."

 A: Well, if I haven't made some calculation error (or overlooked something), it should be like this:
Fix a planar inertial coordinate system $O \, \vec{e}_1 \, \vec{e}_2$ with respect to which the rod is moving. Let $G$ be the center of mass of the rod and let the rod be of length $2\,l$. Attach to the rod a coordinate system $G \, \vec{E}_1 \, \vec{E}_2$ so that the unit vector $\vec{E}_1$ is aligned with the rod, while the unit vector  $\vec{E}_2$ is perpendicular to the rod. Hence, the coordinate system $G \, \vec{E}_1 \, \vec{E}_2$ moves with the rod, relative to the fixed inertial system $O \, \vec{e}_1 \, \vec{e}_2$. At each moment of time $t$ the rod, and hence the vector $\vec{E}_1$ forms (time-changing) angle $\theta = \theta(t)$ with the fixed $O\,\vec{e}$ axis. Consequently,
\begin{align}
&\vec{E}_1 = \vec{E}_1(\theta) = \,\,\,\,\, \cos(\theta)\, \vec{e}_1 + \sin(\theta)\, \vec{e}_2\\
&\vec{E}_2 = \vec{E}_2(\theta) = - \sin(\theta)\, \vec{e}_1 + \cos(\theta)\, \vec{e}_2
\end{align}
Let $\vec{e}_3$ be a vector perpendicular to the plane of $O \, \vec{e}_1 \, \vec{e}_2$. Set
$$\vec{r} \, = \, x \,\vec{e}_1  \, + \, y \, \vec{e}_2$$ to be the position vector pointing from the origin $O$ to the rod's center of mass $G$. Let the moment of inertia of the rod be $$I = \frac{m (2l)^2}{12} = \frac{m l^2}{3}$$
The non-holonomic restriction implies that there is a force acting on the rod, applied to the end, let's say, $l\vec{E}_1$ and always perpendicular to the rod, so $$\vec{f} \, = \, \lambda \, \vec{E}_2(\theta)$$ for some coefficient $\lambda$ which we will determine from the non-holonomic constraint. The end of the rod with the blade has position
$$\vec{r} \, + \, l \, \vec{E}_1(\theta)$$ and its velocity is therefore
$$\vec{v} \, = \, \frac{d \vec{r}}{dt} \, + \, l \frac{d\vec{E}_1}{dt}(\theta)$$
However, the angular velocity of the rod is $$\vec{\omega} \, = \, \frac{d\theta}{dt} \, \vec{e}_3$$ so $$\frac{d\vec{E}_1}{dt}(\theta) \, = \, \frac{d\theta}{dt} \, \vec{e}_3 \times \vec{E}_1(\theta) \, = \, \frac{d\theta}{dt} \, \vec{E}_2(\theta)$$
$$\vec{v} \, = \, \frac{d \vec{r}}{dt} \, + \, l  \, \frac{d\theta}{dt} \, \vec{E}_1(\theta)$$ The non-holonomic constraint imposes the equality $$\vec{v} \cdot \vec{E}_2(\theta) = 0$$
which more explicitly is
$$\left(\frac{d \vec{r}}{dt}\cdot \vec{E}_2(\theta)\right) \, + \, l  \, \frac{d\theta}{dt} \, \left( \, \vec{E}_2(\theta)  \cdot \vec{E}_2(\theta)\, \right) \, = \, 0$$ Since, however, $\vec{E}_2(\theta) \cdot \vec{E}_2(\theta)\, =  \, 1$ we get the holonomic constraint
$$l \, \frac{d\theta}{dt} \, = \, -\, \left(\frac{d \vec{r}}{dt}\cdot \vec{E}_2(\theta)\right) $$
If we combine the equations of motion with this constraint, we get the system
\begin{align}
& m \frac{d^2 \vec{r}}{dt^2} \, = \, \vec{f}\\
&\\
& I \frac{d \vec{\omega}}{dt} \, = \, l\, \vec{E}_1(\theta) \times \vec{f}\\
&\\
&l \, \frac{d\theta}{dt} \, = \, -\, \left(\frac{d \vec{r}}{dt}\cdot \vec{E}_2(\theta)\right) 
\end{align}
and since, as already discussed, $\vec{f} \, = \, \lambda \, \vec{E}_2(\theta) $
\begin{align}
& m \frac{d^2 \vec{r}}{dt^2} \, = \, \lambda \, \vec{E}_2(\theta)\\
&\\
& I \frac{d^2 \vec{\theta}}{dt^2} \, \vec{e}_3 \, = \, \lambda \, l\, \vec{E}_1(\theta) \times \vec{E}_2(\theta) \, = \, \lambda \, l\, \vec{e}_3\\
&\\
&l \, \frac{d\theta}{dt} \, = \, -\, \left(\frac{d \vec{r}}{dt}\cdot \vec{E}_2(\theta)\right) 
\end{align}
which becomes
\begin{align}
& m \frac{d^2 \vec{r}}{dt^2} \, = \, \lambda \, \vec{E}_2(\theta)\\
&\\
& I \frac{d^2 \vec{\theta}}{dt^2} \, = \, l \, \lambda\\
&\\
&l \, \frac{d\theta}{dt} \, = \, -\, \left(\frac{d \vec{r}}{dt}\cdot \vec{E}_2(\theta)\right) 
\end{align}
To have a well-defined system of ODEs we need to determine the unknown $\lambda$. This is done by differentiating the last equation with respect to $t$, then plugging the first two in the differentiated result, and solving for $\lambda$.
$$\frac{d}{dt} \, \Big(\, l \, \frac{d\theta}{dt}\,\Big) \, = \, -\, \frac{d}{dt} \left(\frac{d \vec{r}}{dt}\cdot \vec{E}_2(\theta)\right) $$
\begin{align}
l \, \frac{d^2\theta}{dt^2} \, =& \, -\, \left(\frac{d^2 \vec{r}}{dt^2}\cdot \vec{E}_2(\theta)\right) \, -\, \left(\frac{d \vec{r}}{dt}\cdot \frac{d\vec{E}_2}{dt}(\theta)\right) \, \\
=& \, -\, \left(\frac{d^2 \vec{r}}{dt^2}\cdot \vec{E}_2(\theta)\right) \, -\, \left(\frac{d \vec{r}}{dt}\cdot \frac{d\vec{E}_2}{d\theta}(\theta) \,\right)\frac{d\theta}{dt}\\
=& \, - \, \left(\frac{d^2\vec{r}}{dt^2} \cdot \vec{E}_2(\theta)\right) \, + \, \left(\frac{d\vec{r}}{dt}\cdot \vec{E}_1(\theta) \right) \frac{d\theta}{dt} 
\end{align}
so
\begin{align}
\frac{l^2}{I} \, \lambda \, 
=& \, - \, \frac{1}{m}\left(\vec{E}_2(\theta) \cdot \vec{E}_2(\theta)\right)\, \lambda \, + \, \left(\frac{d\vec{r}}{dt}\cdot \vec{E}_1(\theta) \right) \frac{d\theta}{dt}\\
\frac{3}{m} \, \lambda \, 
=& \, - \, \frac{1}{m}\,\lambda \, + \, \left(\frac{d\vec{r}}{dt}\cdot \vec{E}_1(\theta) \right) \frac{d\theta}{dt}\\
\lambda \, 
=& \, \frac{m}{4}\, \left(\frac{d\vec{r}}{dt}\cdot \vec{E}_1(\theta) \right) \frac{d\theta}{dt}
\end{align}
The complete system of equations of motion is then
\begin{align}
& \frac{d^2 \vec{r}}{dt^2} \, = \, \frac{1}{4}\, \left(\frac{d\vec{r}}{dt}\cdot \vec{E}_1(\theta) \right) \frac{d\theta}{dt} \,\, \vec{E}_2(\theta)\\
&\\
&\frac{d^2 \vec{\theta}}{dt^2} \, = \,  \, \frac{3}{4l}\, \left(\frac{d\vec{r}}{dt}\cdot \vec{E}_1(\theta) \right) \frac{d\theta}{dt}
\end{align}
If you write them componentwise
\begin{align}
&\frac{d^2x}{dt^2} \, =  \, - \, \frac{1}{4} \left(\,\frac{dx}{dt} \cos(\theta) + \frac{dy}{dt} \sin(\theta) \,\right) \frac{d\theta}{dt} \, \sin(\theta)\\
&\frac{d^2y}{dt^2} \, =  \,\,\,\,\,\,\,    \frac{1}{4} \left(\,\frac{dx}{dt} \cos(\theta) + \frac{dy}{dt} \sin(\theta) \,\right) \frac{d\theta}{dt} \, \cos(\theta)\\
&\frac{d^2\theta}{dt^2} \, =  \,\,\,\,\,\,\,   \frac{3}{4l} \left(\,\frac{dx}{dt} \cos(\theta) + \frac{dy}{dt} \sin(\theta) \,\right) \frac{d\theta}{dt}
\end{align}
