4. Lagrangians and Hamiltonians

Legal Notice, Disclaimer, and Book Ordering Information:

Null

Setup

\(Clear["Global`*`"]; 
Off[General::spell]; 
Off[General::spell1]; 
    \)

<< "D:\Mma\Inits\ch4init.m"

4.2 Problems

4.2.1  Lagrangian Problems

Problem 3:  Bead on a Rotating Hoop

Clear["Global`*`"];

Part (a)   Construct the Lagrangian and derive the equations of motion.

x2rRule = 
      Thread[{x, y, z} 
          -> {r Sin[θ[#1]] Cos[ω #1]&, 
            r Sin[θ[#1]] Sin[ω #1]&, 
            r Cos[θ[#1]]&}];

T=(1/2) m (((x'[t])^2)+((y'[t])^2)+((z'[t])^2));V = m g z[t];

L = 
      ExpandAll[
        TrigReduce[Simplify[T - V //. x2rRule]] //. {
            Cos[x_]^2 -> 1 - Sin[x]^2}]

Lag[θ, L]

θSol=ExpandAll[(Flatten //. {g -> ω(c^2) r})[Solve[Lag[θ, L],θ''[t]]]]

energy=ExpandAll[θ'[t] (∂_(θ'[t]))L-L//.{g -> ω(c^2) r}//.{Cos[x_]^2 -> 1 - Sin[x]^2}]

{pRule, hamiltonian, eqMotion} = 
      Hamilton[L, {θ}, {pθ}];

energy == hamiltonian //. {g -> ω(c^2) r} //. pRule

Expand[
      (∂_t) energy //. θSol //. {
          g -> ω(c^2) r}]

Part (b)   Construct the conserved energy and verify that it is conserved.

Veff=energy//.{θ'[t]->0}

dVeff = Factor[(∂_(θ[t]))Veff]

Solve[(dVeff/Sin)[θ[t]] == 0, θ[t]]
        [[1]]

eqRule = 
      θ[t] 
        -> {0, π, ArcCos[-ω((c^2)/(ω^2))], 
          ArcCos[((+ω(c^2))/(ω^2))] + π};

Thread[dVeff == 0 /. eqRule]

Part (c)   Define a potential for the system and find the equilibrium points.  Show that the characteristics of the potential change at ω=ωc, where ωc^2 =g/r.  Graph the potential for the two characteristic types of potentials.

values1 = {r -> 1, ωc -> 1, m -> 1, 
        θ[t] -> θ};

Plot3D[
      Evaluate[Veff //. values1], {θ, -π, π}, {ω, 
        1, 3}, AxesLabel 
        -> {"θ", "ω", "Potential"}];

p1 = 
      Plot[Evaluate[Veff //. values1 /. {ω -> 0.2}], {
          θ, -π, π}];

p2 = 
      Plot[Evaluate[Veff //. values1 /. {ω -> 2.0}], {
          θ, -π, π}];

Show[p1, p2];

Part (d)   Numerically solve for θ[t] and θ'[t] and plot the data.

eq1 = 
      θSol //. {g -> ω(c^2) r} /. 
        Rule -> Equal

eq2=Flatten[{eq1/.{θ''[t]->(dθ)'[t]},θ'[t]==dθ[t]}]

values2 = {ω -> 2, ωc -> 1};

initial = {θ[0] == π, dθ[0] == 0.01};

ndsol1 = 
      NDSolve[
          Flatten[{eq2 //. values2, initial}], {θ[t], dθ[t]}, {
            t, 0, 30}][[1]]

\(θ1[t_] = 
      θ[t] /. ndsol1[[1]]; 

    dθ1[t_] = 
      dθ[t] /. ndsol1[[2]]; \)

Plot[θ1[t], {t, 0, 30}];

Plot[dθ1[t], {t, 0, 30}];

Part (e)   Animate the numerical results.

bead[t_] := 
      Graphics[{PointSize[0.03], 
          Point[{Sin[θ1[t]], Cos[θ1[t]]}]}];

\(frames = 20; 

    Do[Show[bead[t], 
        Graphics[{Hue[ .08], Circle[{0, 0}, 1]}, 
          AspectRatio -> Automatic], 
        PlotRange -> {{\(-2\), 2}, {\(-2\), 2}}], {t, 0, 18, 
        18\/\(frames - 1\)}]\)