Derivation of Lagrange Equations from Newton's Second Law for a Non-holonomic System of Particles I am interested to write down a derivation of Lagrange equations from Newton's second law for a non-holonomic system of particles. Here, I mention my derivation where I am stuck right at the last step.
Consider a system of $N$ particles where their position vectors are written as
$$\mathbf{r}_i=\mathscr{R}_i(q_1(t),\dots,q_M(t),t),\quad i=1,\dots,N\,,\tag{1}$$
where the functions $q_i:\mathbb{R}\to\mathbb{R}$ are called the generalized coordinates which are subjected to holonomic and non-holonomic constraints as below
\begin{align*} 
f_i(q_1(t),\dots,q_M(t),t)&=0,\quad i=1,\dots,C_h\,, \\
g_i(q_1(t),\dots,q_M(t),\dot q_1(t),\dots,\dot q_M(t),t)&=0,\quad i=1,\dots,C_n\,,
\tag{2}
\end{align*}
where $C_h$ and $C_n$ are the number of holonomic and non-holonomic constraints, respectively. Also, if the degree of freedom of the system is $n$ then $n=M-C\ge1$ where $C=C_n+C_h$ is the total number of constraints. Using the chain rule of differentiation we have
\begin{align*}
\mathscr{\dot R}_i := \mathbf{v}_i &= \mathbf{v}^*_i+\frac{\partial\mathscr{R}_i}{\partial t},\quad i=1,\dots,N\,, \\
\mathbf{v}^*_i&:=\sum_{j=1}^{M}\frac{\partial \mathscr{R}_i}{\partial q_j}\dot q_j\,,\tag{3}
\end{align*}
where we defined the virtual velocity of a particle by $\mathbf{v}^*_i$. Also, from Newton's second law we have
$$\mathbf{F}_i=m \mathbf{a}_i\tag{4},\quad i=1,\dots,N\,.$$
Multiplying both sides of $(4)$ by $\mathbf{v}^*_i$, summing over the number of particles $N$ and interchanging the the order of summations we get
$$\sum_{j=1}^{M}\sum_{i=1}^{N}(\mathbf{F}_i-m\mathbf{a}_i)\cdot\frac{\partial \mathscr{R}_i}{\partial q_j}\dot q_j=0\,.\tag{5}$$ 
Then using the following definitions and identities
\begin{align*}
Q_j&:=\sum_{i=1}^{N}\mathbf{F}_i\cdot\frac{\partial \mathscr{R}_i}{\partial q_j},\quad j=1,\dots,M\,, \\
S_j&:=\sum_{i=1}^{N}m\mathbf{a}_i\cdot\frac{\partial \mathscr{R}_i}{\partial q_j}=\frac{d}{dt}\frac{\partial T}{\partial \dot q_j}-\frac{\partial T}{\partial q_j},\quad j=1,\dots,M\,, \\
T&:=\sum_{i=1}^{N}\frac{1}{2}m\mathbf{v}_i\cdot\mathbf{v}_i\,,
\tag{6}
\end{align*}
Eq. $(5)$ reduces to
$$\sum_{j=1}^{M}(Q_j-S_j)\dot q_j=0.\tag{7}$$
If there were no constraint equations at all, either holonomic or non-holonomic as mentioned in Eq.$(2)$, then the functions $q_i$ were linearly independent and from this we could conclude that the functions $\dot q_i$ are also linearly independent. Then Eq.$(7)$ would result in the well known form of Lagrange equations $S_j=Q_j$. But here is my question, what if there are constraint equations like Eq.$(2)$. Note that sometimes we are inclined not to eliminate the holomonic constraints by using a transformation. So I am insisting to have both holonomic and non-holonomic constraints at the same time.

As the functions $\dot q_i$ are not (linearly) independent in this case, I am wondering that how the last step works here?

If the non-holonomic constraints are linear in terms of generalized velocities
\begin{align*}
&g_i(q_1(t),\dots,q_M(t),\dot q_1(t),\dots,\dot q_M(t),t)=\\
&\sum_{j=1}^{M}a_{ij}(q_1(t),\dots,q_M(t),t)\dot q_j(t)+b_i(q_1(t),\dots,q_M(t),t)=0,\quad
\tag{8}
\end{align*}
then we call them quasi non-holonomic. In this case, I know that the final result should be
$$\frac{d}{dt}\frac{\partial T}{\partial \dot q_j}-\frac{\partial T}{\partial q_j}=Q_j+\sum_{i=1}^{C_h}\lambda_i\frac{\partial f_i}{\partial q_j}+\sum_{i=1}^{C_n}\mu_i\frac{\partial g_i}{\partial \dot q_j},\quad j=1,\dots,M\,,\tag{9}$$
where $\lambda_i$ and $\mu_i$ are some functions of time which are called Lagrange multipliers.
A simple observation is that every holonomic constraint can be written in the form of a quasi non-holonomic constraint, that is
\begin{align*}
&\\
&\sum_{j=1}^{M}\frac{\partial f_i}{\partial q_j}(q_1(t),\dots,q_M(t),t)\dot q_j(t)+\frac{\partial f_i}{\partial t}(q_1(t),\dots,q_M(t),t)=0.\quad
\tag{10}
\end{align*}
At the first step, it seems reasonable to establish an argument when all of the constraints are quasi non-holonomic. I have posted a related mathematical question on Mathematics SE in this regard.  Interested reader can take a look at it.
 A: The full question (v13) is rather broad, but here are some comments/feedback:


*

*Traditionally virtual displacements are frozen in time $t$. The $t$-differentiations in eqs. (3) & (5) are misleading at best (depending on what the notation $t$ is supposed to represent).

*Eq. (4) is Newton's 2nd law if ${\bf F}_i$ denotes the total force on the $i$th point particle, i.e. a sum of "applied" and constraint forces. One of the main points is to try to eliminate the constraint forces from the formalism, at least for the holonomic constraints. OP seems to have made no progress in this.

*It should be noted that a $\partial g / \partial \dot q_j$-term in OP's final eq. (8) is not appropriate for a general non-holonomic constraint $g$, but only for a so-called semi-holonomic constraint, 
$$ g(q,\dot{q},t)~\equiv~\sum_j a_j(q,t)\dot{q}^j+a_t(q,t)~\approx~0, $$
which by definition is an affine function of $\dot{q}^j$.

*In case OP is following Ref. 1, note that the treatment of Lagrange equations for non-holonomic constraints in Ref. 1 is inconsistent with Newton's laws, and has been retracted on the errata homepage for Ref. 1. See Ref. 2 for details. 

*In case OP is following Ref. 3, note that Ref. 3 only deals with holonomic constraints in Chapter 2. (This is explicitly mentioned on the middle of p. 50.) Ref. 3 tentatively introduces non-holonomic constraints in the beginning of Section 3.1.2, only to later rejects the approach as unphysical.
References:


*

*H. Goldstein, Classical Mechanics; 3rd ed; Section 2.4. Errata homepage. (Note that this criticism only concerns the treatment in the 3rd edition; the results in the 2nd edition are correct.)

*M.R. Flannery, The enigma of nonholonomic constraints,
Am. J. Phys. 73 (2005) 265.

*J.V. Jose & E.J. Saletan, Classical Dynamics: A Contemporary Approach, 1998; Sections 2.1.1, 2.2.1 & 3.1.2.
A: The problem is that you can't rearrange the formulae defining momenta to express the $\dot{q}_i$ in terms of them and the $q_j$, because of phase space constraints of the form $f(q,\,p)=0$. The trick is to identify these functions and add multiples of them to the Hamiltonian, viz. https://en.wikipedia.org/wiki/Dirac_bracket
A: Backwards, and with boring rectilinear coordinates:
$$ \begin{eqnarray}\frac{\partial \mathcal L}{\partial x_i} &=& \frac{d}{dt} \frac{\partial \mathcal L}{\partial \dot x_i} \\
\mathcal L = T - U &=& \frac{m}{2} \sum_i \dot x_i^2-  \int_C \mathbf F \cdot d \mathbf x \\
\frac{\partial \mathcal L}{d x_i} = - \frac{\partial U}{\partial x_i} = - \frac{\partial}{\partial x_i} \int_C \mathbf F \cdot d\mathbf x &=& \frac{d}{dt} \frac{\partial \mathcal L}{\partial \dot x_i} =  \frac{d}{dt} \frac{\partial T}{\partial \dot x_i} =  \frac{d}{dt} m \dot x_i  \\
\mathbf F &=& m \mathbf{\ddot x}
\end{eqnarray}$$
If you're so inclined, the forwards derivation flows the opposite direction and the least intuitive steps are from $m \mathbf{\ddot x}$ back to the derivatives of $T$ term, but all the steps are valid. 
Generalized coordinates are formally equivalent to rectilinear coordinates even if quirky terms like effective potentials arise (when the kinetic energy has derivatives with respect to position in addition to velocities). If you want to be explicit about it, Wikipedia does a good job.
It sounds like you're also looking for information about non-conservative systems. Please examine this link. It also sounds like you're curious about the constraints, see here.
A: After a while, I noticed that I have thought so badly due to bad books that I have read and bad explanation that I have been given.
Let me start by deriving the Lagrange equations for the motion of a single particle. First, choose a set of generalized coordinates. For example, consider a particle moving on a spherical surface whose radius is varying with time. You can choose the generalized coordinates to be just the usual Cartesian coordinates $(x,y,z)$ subjected to the holonomic constraint
$$f(x(t),y(t),z(t),t)=x(t)^2+y(t)^2+z(t)^2-R^2(t)=0,$$
where $R(t)$ is the radius of the sphere. Another way is to choose the spherical coordinates $(r,\theta,\phi)$ as the generalized coordinates
\begin{align*}
x&=r\sin\theta\cos\phi \\
y&=r\sin\theta\sin\phi \\
z&=r\cos\theta
\tag{1}
\end{align*}
subjected to the constraint $r=R(t)$, or we may just choose $(\theta,\phi)$ as the generalized coordinates with $r$ being replaced by $R(t)$ in $(1)$. You can also choose cylindrical coordinates $(\rho,\phi,z)$ as the generalized coordinates subjected to the constraint $\rho^2+z^2=R^2(t)$. These were just simple examples of what the equation
$$\mathbf{r}=\mathscr{R}(q_1(t),\dots,q_M(t),t), \tag{2}$$
means. It is wise to pick that number of generalized coordinates which do not exceed the maximum degree of freedom that the particle can have so $M\le3$. Note that we do not care about the fact that the system is constrained or not when we choose our generalized coordinates. Only after choosing the generalized coordinates, we seek what the governing constraints are. It is natural to choose those generalized coordinates, which make our description of motion and the solution of the problem easy. That's all.
Now, write down Newton's second law for the particle and multiply both sides with the vector $\frac{\partial \mathscr{R}}{\partial q_j}$ to obtain
$$(\mathbf{F}-m\mathbf{a})\cdot \frac{\partial \mathscr{R}}{\partial q_j}=0\,.\tag{3}$$
Using the definitions and identities mentioned in Eq.$(6)$ of the above question this simply turns out to be
$$Q_j-S_j=0.\tag{4}$$
But, there is still an important question. How many linearly independent equations do we get!? A close look at $(3)$ reveals that this depends on how many linearly independent vectors $\frac{\partial \mathscr{R}}{\partial q_j}$ we have. This serves as a motivation to require that the generalized coordinated be chosen such that all of the vectors $\frac{\partial \mathscr{R}}{\partial q_j}$ be linearly independent. In the aforementioned example of the moving particle over the sphere, if we choose spherical coordinates, then these vectors are simply parallel to basis vectors of spherical coordinates, which are indeed linearly independent. Keeping this in mind, we shall get $M$ linearly independent equations.
Note that our discussion never depends on whether the motion is constrained or not. It is just another story. The corresponding forces of kinematical constraints show them-selves in $Q_j$. For the case of holonomic constraints, one can show that if a set of generalized coordinates is chosen such that the holonomic constraint is satisfied identically, then the contribution of its corresponding forces to $Q_j$ vanishes. This corresponds to our choice $(\theta,\phi)$ in the above example. In conclusion, depending on the choice of generalized coordinates, the corresponding constraint forces may or may not enter the Lagrange equations.
With slight modifications of the procedure in the above question, you can derive the Lagrange's equations for a system of particles. At the first step, you get
$$\sum_{i=1}^{N}(\mathbf{F}_i-m\mathbf{a}_i)\cdot\frac{\partial \mathscr{R}_i}{\partial q_j}=0,$$
which is again equivalent to $Q_j-S_j=0$. Note that no summation on $j$ is really needed! The tricky point is to understand how many linearly independent equations you are left with, which I leave as an exercise for you.
Note that this method tells us that Lagrange's equations are just well-chosen linear combinations of the original equations of motion obtained by Newton's second law and nothing more! It is well-chosen as it eliminates some forces which we are not interested to know about.
Here is a nice reference which I strongly recommend for interested readers.
Classical Dynamics: A Contemporary Approach By Jorge V. José and‎ Eugene J. Saletan.
Furthermore, for a more solid background on this stuff, you may want to study differential geometry. Specifically, a thorough understanding of the Krush-Kahn-Tucker theorem is the key to understanding constrained dynamics.

Disclaimer. These were my understanding of the subject and it may contain some mistakes. So, I would appreciate any helpful and detailed criticism.
