Equation of motion through the Lagrangian with Lagrange multipliers I ask for advice, cause I'm a little confused.
We have such a Lagrangian:
$L=\frac{1}{2}m(\dot{x}^2+\dot{y}^2)-\lambda(x+xy+y-1)$
Here $\lambda(x+xy+y-1)$ is the constraint on the phase variables.
I need to derive the equation of motion given the constraints and solve them numerically with the help of NDSolve.
We do this in accordance with the classic formula:
$\frac{d}{dt}(\frac{dL}{d\dot{q}})-\frac{dL}{dq}=0$
Where $q=[x,y]$ are generalized coordinates.  I'm not sure about the Lagrange multiplier as a generalized coordinate.
Clear["Derivative"]

ClearAll["Global`*"]

T = 1/2 m (x'[t]^2 + y'[t]^2);(*Kinetic Energy*)

f = \[Lambda] (x[t] + x[t] y[t] +y[t] - 1);(*Constraint*)

L = T - f;(*Lagrangian*)

D[D[L, x'[t]], t] - D[L, x[t]];

D[D[L, y'[t]], t] - D[L, y[t]];

D[D[L, \[Lambda]'[t]], t] - D[L, \[Lambda][t]];

Question: how are the Lagrange multipliers included in this system when compiling the ODE system and numerically solving it?
Maybe this help?
https://farside.ph.utexas.edu/teaching/336k/lectures/node90.html
https://www.sciencedirect.com/science/article/abs/pii/0045782588900850
 A: If the underlying Lagrangian is just kinetic energy (i.e. there is no potential energy in the system) and you need to stay on the curve $F=0$, then the Newton equations for the constrained Lagrangian dynamics become just $m\ddot{x}-\lambda(t) F_x=0,m\ddot{y}-\lambda(t) F_y=0$, so the force points parallel to $\nabla F$ but with a coefficient that depends on time, and the requirement that the constraint is satisfied all the time determines the evolution of $\lambda$. In order for these equations to be consistent, the initial velocity will need to be tangent to the curve, but in this case you can solve the two differential equations together with the algebraic equation to describe the evolution.
Note that this constrained Lagrangian dynamics physically corresponds to "making up the right force that keeps you on the specified curve" (so it is a bit less magical than it looks at first glance). For a classic example you can think of moving on the circle $x^2+y^2=r^2$ starting at $(r,0)$ with initial velocity $(0,v)$. From this you can recover the familiar $\mathbf{a}=-\frac{v^2}{r} \mathbf{e}_r$.
A perhaps more physical view of a situation like this in the absence of any meaningful potential energy is to just take some strictly convex function $G$ minimized at zero (e.g. $G(x)=x^2$). and run the evolution throughout with $m(\ddot{x},\ddot{y})^T=-\nabla G(F(x,y))$.
A: This is a DAE system of index 2 or 3. You obtain an ODE system by adding derivatives of the equation, mainly the constraint, to the system.
$$
x+xy+y=1,\\
\dot x(1+y)+\dot y(1+x)=0,\\
\ddot x(1+y)+\ddot y(1+x)+2\dot x\dot y=0
$$
From the last equation eliminate the second derivatives to get an equation for the Lagrange multiplier. With another derivative (thus index 3) you get an ODE for the multiplier.
$$
m\ddot x=-λ(1+y)\\
m\ddot y=-λ(1+y)\\
\implies 
-λ[(1+x)^2+(1+y)^2]+2m\dot x\dot y=0\\
λ=\frac{2m\dot x\dot y}{(1+x)^2+(1+y)^2}
$$
