0
$\begingroup$

I ask for advice, I'm a little confused. I have such a Lagrangian.

$L=\frac{1}{2}m(\dot{x}^2+\dot{y}^2)-\lambda(2x^2+3y^2-1^2)$

Here $\lambda(2x^2+3y^2-1^2)$ 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] (2 x[t]^2 + 3 y[t]^2 - 1^2);(*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?

EDIT (with correct initial conditions and another constraint):

Clear["Derivative"]

ClearAll["Global`*"]

q = Function[t, {x[t], y[t]}];(*coordinates*)

zwang = x[t] + x[t] y[t] + y[t] - 1;(*constraints*)

Solve[zwang == 0, y[t]];

Q = \[Lambda][t] D[zwang, {q[t]}];(*gen.forces*)

L = 1/2 (x'[t]^2 + y'[t]^2);

eqn = D[D[L, {q'[t]}], t] - 
   D[L, {q[t]}] - \[Lambda][t] D[zwang, {q[t]}];


sol = Solve[
    Join[eqn, {D[zwang, {t, 2}]}] == 0, {x''[t], 
     y''[t], \[Lambda][t]}][[1]];

x0 = 1/2;

y0 = 1/3;

s = NDSolve[{sol[[1, 1]] == sol[[1, 2]], sol[[2, 1]] == sol[[2, 2]], 
    sol[[3, 1]] == sol[[3, 2]], x[0] == x0, y[0] == y0, x'[0] == 0, 
    y'[0] == 0}, {x, y, \[Lambda]}, {t, 500}];

Plot[Evaluate[zwang /. s], {t, 0, 500}, PlotRange -> All]
$\endgroup$
11
  • 1
    $\begingroup$ [Lambda] is a constant, the derivative would be zero. Instead you need to use the equation of the constrain: x^2+y^2-1==0 $\endgroup$ Oct 14 '21 at 8:14
  • $\begingroup$ @DanielHuber Yes, I agree. As for the second remark, constraint equation does not contain derivatives. This means that the system of ODE will include algebraic constraints, with such systems, as far as I know, acceptable in Mathematica. Is that correct? $\endgroup$
    – dtn
    Oct 14 '21 at 8:17
  • $\begingroup$ Which of the constraints is the correct one: 2 x[t]^2 + 3 y[t]^2 ==1^2 or x[t]^2 + y[t]^2 ==1^2 ? $\endgroup$ Oct 14 '21 at 8:24
  • 1
    $\begingroup$ Have a look at: farside.ph.utexas.edu/teaching/336k/lectures/node90.html $\endgroup$ Oct 14 '21 at 9:03
  • 1
    $\begingroup$ @DanielHuber Thanks for this very interesting link! $\endgroup$ Oct 14 '21 at 12:35
2
$\begingroup$

As far as I remember Lagrangeformalism (classic form) is only valid for generalized coordinates, which fullfill the constraints . In your example x[t],y[t] aren't generalized coordinates!

For the constraint 2 x[t]^2 + 3 y[t]^2 ==1^2 possible generalized coordinate q[t] would be x[t]->Cos[q[t]]/Sqrt[2],y[t]->Sin[q[t]]/Sqrt[3]

2 x[t]^2 + 3 y[t]^2 ==1^2 /.{x[t]->Cos[q[t]]/Sqrt[2],y[t]->Sin[q[t]]/Sqrt[3]}// Simplify
(*True*)

L= 1/2 m (3 x'[t]^2 + y'[t]^2) /. {x -> ( Cos[q[#]]/Sqrt[2] &), y  -> (Sin[q[#]]/Sqrt[3] &)} //Simplify
(*-(1/24) m (-11 + 7 Cos[2 q[t]]) Derivative[1][q][t]^2*)

The equations of motion follow to

D[D[L, q'[t]], t] - D[L, q[t]]==0
(*7/12 m Sin[2 q[t]] Derivative[1][q][t]^2 -1/12 m (-11 + 7 Cos[2 q[t]]) (q^\[Prime]\[Prime])[t]*)

Hope it helps!

addendum Lagrange with constraints (thanks @DanielHuber's comment )

q = Function[t, {x[t], \[CurlyPhi][t]}]; (*coordinates*)
zwang = x[t] - r \[CurlyPhi][t] (* constraints *)
Q = \[Lambda][t] D[zwang, {q[t]}] (* gen. forces*)

L = 1/2 m x'[t]^2 + 1/2 \[Theta] \[CurlyPhi]'[t]^2 +m g Sin[\[Alpha]]  x[t]
eqn = D[D[L, {q'[t]}], t] -D[L, {q[t]}] - \[Lambda][t] D[zwang, {q[t]}]

From equations eqn and zwang the Lagrange multiplier \[Lambda][t]is evaluated, substitution gives the equations of motion

sol = Solve[Join[eqn, {D[zwang, {t, 2}]}] ==0, {x''[t], \[CurlyPhi]''[t], \[Lambda][t]}][[1]]
(*{(x^\[Prime]\[Prime])[t] -> (g m r^2 Sin[\[Alpha]])/(m r^2 + \[Theta]), (\[CurlyPhi]^\[Prime]\[Prime])[t] -> (g m r Sin[\[Alpha]])/(m r^2 + \[Theta]), \[Lambda][t] -> -((g m \[Theta] Sin[\[Alpha]])/(m r^2 + \[Theta]))}*)

Q /. sol (* Reibkraft und -moment*)

 
$\endgroup$
10
  • 2
    $\begingroup$ I answered the question you asked for. In this new example substitute x[t]-> (1 - y[t]^2)/(1 + 3 y[t])) , Lagrangian follows to L=(m (5 + 12 y[t] + 38 y[t]^2 + 60 y[t]^3 + 45 y[t]^4) Derivative[1][y][ t]^2)/(1 + 3 y[t])^4 with generalized coordinate y[t] $\endgroup$ Oct 14 '21 at 8:53
  • 1
    $\begingroup$ I tried to solve the constraints Solve[x + 3 x y + y^2 \[Minus] 1 == 0, x]. $\endgroup$ Oct 14 '21 at 8:55
  • 1
    $\begingroup$ My modified answer shows how to proceed! You could check your results on your own I think. $\endgroup$ Oct 16 '21 at 11:04
  • 1
    $\begingroup$ In your edit the initial conditions doesn't fullfill the constraint: zwang /. t -> 0 /. {x[0] -> 1/2, y[0] -> 1/4} (*-1/8 !=0 *) $\endgroup$ Oct 18 '21 at 8:44
  • 1
    $\begingroup$ Fine! But still in your edit you are using wrong initial conditions... $\endgroup$ Oct 18 '21 at 9:20

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.