Legal Notice, Disclaimer, and Book Ordering Information:

(* Legal Notice, Disclaimer, and Book Ordering Information:

Last Modified: FIO June 8, 2002

Mathematica For Physics: 2nd Edition

Robert L. Zimmerman and Fredrick I. Olness

WebSite:  
  www.physics.smu.edu/~olness
  darkwing.uoregon.edu/~phys600/

MathSource Number: 0206-862

ISBN 0-201-53796-6
For ordering information, call 1-800-282-0693
Addison-Wesley Publishing Company

Communication with the authors:
Fredrick I. Olness:  olness@mail.physics.smu.edu
Robert L. Zimmerman: bob@zim.uoregon.edu

Copyright 2002, Addison-Wesley Publishing Company, Inc. The material in
this file may be distributed freely  so long as the content remains
unchanged, and the copyright and  reference notices are included. The
publisher grants permission for the noncommercial use of these programs
and program segments.  All other uses require the prior written consent
of the publisher.    

All rights are reserved with regard to the  material in the
"Mathematica For Physics" text, and may not be reproduced, stored in a
retrieval system, or transmitted,  in any form or by any means,
electronic, mechanical, photocopying,  recording, or otherwise, without
the prior written permission  of the publisher. 

The programs and applications presented in this book have been included
for  their instructional value. They have been tested with care but are
not  guaranteed for any particular purpose. The publisher does not
offer any  warranties or representations, nor does it accept any
liabilities with  respect to the programs or applications. 

 *)

Contains correction on:  
    2nd Edition:Chpt.6,Problem 6.2.6 p.334:Spring Pendulm. 
    2nd Edition: Chpt.6, Problem 6.3.3 p.174: Sphere Rolling on a Fixed Sphere.
(FIO, 17 Jan 2003)  

Contains correction on:  
    2nd Edition:Chpt.6,Problem 6.4.2 p.351: Nonlinear Oscillator  (FIO, 19 Apr 2004)  
    
Contains correction on:  
    2nd Edition:Chpt.6,Problem 6.2.4 p.319:  Springs Mounted on Top of a Carriage  (FIO, 20 Sept. 2006)

6. Lagrangians and Hamiltonians

6.1 Introduction

Off[General::spell];  
Off[General::spell1];

6.2 Lagrangian Problems without Lagrange Multipliers

Overview of  Lagrangian Problems without Lagrange Multipliers

Initialization for Lagrangian Problems

Needs["Calculus`VariationalMethods`"]

?EulerEquations

L = m/2  x'[t]^2  - k/2  x[t]^2

EulerEquations[L, x[t], t] // ExpandAll

?FirstIntegrals

FirstIntegrals[L, x[t], t] // ExpandAll

Problem 1: Particle Sliding on a Movable Incline

Solution

Clear["Global`*"]

Part a

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

change = {x -> ((x1[#] + (L0 - s[#]) 
                      Cos[θ]) &),    
             y -> (((L0 - 
                        s[#]) Sin[θ]) &)};

L=T-V /.change //Simplify

eq1 = EulerEquations[L, {x1[t], s[t]}, t]

eqMotion=Solve[eq1, {x1''[t],s''[t]}] //Flatten

Part b

xmomentum = 
      M x1'[t] + m x'[t] /. change  // Simplify

D[xmomentum, t]  /. eqMotion // Simplify

energy = T + V /. change // Simplify

D[energy, t] /. change /. eqMotion // Simplify

FirstIntegrals[L, {x1[t], s[t]}, t]

Part c

Sol = DSolve[
{( 
                  eqMotion /. Rule -> Equal), 
s[
                    0] == 0, s'[0] == 0, 
                x1[0] == 0, x1'[0] == 0} // Flatten, {x1[t], s[t]}, 
            t  ] // Simplify  // Flatten

(((s[t]  /. Sol)   == L0  // 
                  Simplify
 // 
                Solve[# /. t^2 -> x, x] &
) /. 
            x -> t^2  // Flatten   // FullSimplify )

Problem 2: Bead Sliding on a Rotating Wire

Remarks and outline

Solution

Clear["Global`*"]

Part a

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

ψRule = θ ->  (π - ψ) ;

T=(1/2) m (((x'[t])^2)+((y'[t])^2)+((z'[t])^2));
V = 
        m g z[t];
L = T - V  /. 
            x2rRule  /. ψRule  // Simplify

eq1 = EulerEquations[L, r[t], t]

rSol = ( Solve[eq1, r''[t]][[1, 1]] // ExpandAll ) //. 
        Rule -> Equal

eq2 =  DSolve[ {rSol, 
                  r[0] == 0, r'[0] == v0} , r[t], t][[1, 
                1]]  // Simplify  // PowerExpand   // 
        ExpandAll

rRule= MapAt[ Collect[#, {Exp[-x_]  }]  &,  
        eq2,  2]

Part b

tRule = ( Solve[Len == r[t] /. rRule , t , 
                InverseFunctions -> True][[2]]  // 
            Simplify
 // 
          PowerExpand 
)

Series[t /. tRule, {ψ, 0, 0} ]  // 
        Normal//Simplify

Part c

values = {ψ -> π/13, g -> 1, ω -> 2, 
          v0 -> 0};

position = {x[t], y[t], z[t]} //. x2rRule //. 
              rRule /. ψRule //. values;

p1 = 
        ParametricPlot3D[
          position // Evaluate, {t, 0, 6}, 
Boxed -> 
            True, 
Axes -> 
            False, 
BoxRatios -> {1, 1, 
              1}, 
DisplayFunction -> 
            Identity];

discretePoints = 
        Table[position, {t, 0, 6, 0.5}];
discreteGraphics = 
        Graphics3D[{PointSize[0.03], Point /@ discretePoints}];

Show[discreteGraphics, p1, BoxRatios -> {1, 1, 1}];

Problem 3: Bead on a Rotating Hoop

Remarks and outline

Solution

Clear["Global`*"]

Part a

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

ψRule = {θ -> (π -  ψ[#] & )};

T=(1/2) m (((x'[t])^2)+((y'[t])^2)+((z'[t])^2))//.x2rRule /.ψRule//Simplify

V = m g z[t] //. x2rRule  /. ψRule // Simplify

L=(T - V)

eq1 = EulerEquations[L,  ψ[t], t]

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

energy = (ψ'[t] D[L, ψ'[t]] - 
              L) /. {g -> ω(c^2) r} // Simplify

D[energy, t] /. ψSol /. {g -> ωc^2 r}   // 
      Simplify

FirstIntegrals[L,  ψ[t], 
          t] /. {g -> ω(c^2) r} // Simplify

Part b

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

dVeff= D[Veff,ψ[t]]

sol=Solve[dVeff  == 0, ψ[t], 
        InverseFunctions -> True ]

eqRule=Join[{{ψ[t] -> π}}, 
        sol]

dVeff == 0 //. eqRule

values1 = {r -> 1,  ωc -> 1,  
          m -> 1};

Plot3D[
        Veff /. values1 // 
          Evaluate, 
{ψ[
            t], -π, π}, {ω, 1, 
          3}, 
ColorFunction -> 
          Hue , 
Ticks -> {π Range[-1, 1, 
                1/2], Automatic,  Automatic }];

Plot[(Veff //. 
                    values1 /. ω -> #) &   /@   {2, 
              0.2} // Evaluate, 
{ψ[
            t], -π, π}, 
Ticks -> {{  -π  , 0,  π}, 
            Automatic}, 
Epilog -> {
Hue[ .6], AbsolutePointSize[6], 
            Map[Point,  { ψ[t] , Veff} /. eqRule  /. 
                  values1 /. ω -> 
                  2] }
];

Part c

\(\(pt1 = \(ContourPlot[\ \(\(energy /. values1\) /. {m -> 2, 
                      r -> 1\ }\) /. \((ω -> #)\)\ \  // 
                Evaluate, \ {ψ[t], \(-4\) π, \ 
                4\ \ π}, {\(ψ'\)[t], \(-3\), 
                3}\ , \ 
FrameTicks -> {π\ Range[\(-4\), 4, 2], Automatic}, \ 
ColorFunction -> 
                Hue, \ \ 
PlotPoints -> 
                100, \ 
FrameLabel -> {"ψ", "ψ'"}\ , 
DisplayFunction -> 
                Identity\ ] &\)\ \  /@ \ \ {2, 0.2\ };\)\)

Show[GraphicsArray[ pt1 ]];

cylind =  
        ParametricPlot3D[{v, Sin[u], 
            Cos[u]}, 
{u, -π, π}, {v, -3, 
            3}, 
DisplayFunction -> 
            Identity ];

Show[ 
        cylind, 
Graphics3D[
{AbsoluteThickness[
              3], 
Cases[Graphics[pt1[[1]]], _Line, 
                Infinity] /. 
{x_ , y_} :> {y, 
                  1.02  Cos[x ], 
                  1.02  Sin[
                      x] }}
] , 
ViewPoint -> {1, 1, 
            1}, 
DisplayFunction -> $DisplayFunction];

Part d

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

eqRule /. values2 // N

{g → ωc2 r}  // FullForm

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

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

ndsol1 = NDSolve[{eq2 //. values2, ψ[0] == 1.1318, 
                  dψ[0] == #} // Flatten, {ψ[t], 
                dψ[t]}, {t, 0, 10}][[1]] &   /@   {  1, 
          2.5}

ψ[1][t_]= ψ[t] /.ndsol1[[1,1]];
dψ[1][t_]= dψ[t]/.ndsol1[[1,2]];
 ψ[2][t_]= ψ[t] /.ndsol1[[2,1]];
dψ[2][t_]= dψ[t]/.ndsol1[[2,2]];

pt2 = Plot[  {dψ[#][t], ψ[#][t]}  // 
                Evaluate, {t, 0, 
                10}, 
PlotStyle -> {Dashing[{0.04} ], 
                  AbsoluteThickness[
                    2] }, 
Ticks -> {  Automatic, π/2   Range[1, 12, 3] },  
DisplayFunction -> 
                Identity  ]  &  /@   Range[ 2] ;

Show[GraphicsArray[pt2]];

pt3 = ParametricPlot[  { ψ[#][t],  dψ[#][
                    t]}  // Evaluate, {t, 0, 
                10}, 
PlotStyle ->  
                AbsoluteThickness[
                  2],  
DisplayFunction -> 
                Identity ] &   /@   {1, 2}   ;

Show[ pt1[[1]], pt3,  
        DisplayFunction -> $DisplayFunction];

Part e

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

frames = 2;
\(Do[
Show[{bead[t], 
          Graphics[{Hue[0.1], Circle[{0, 0}, 1]}, 
            PlotRange → {{\(-2\), 2}, {\(-2\), 
                  2}}, 
AspectRatio -> 
              1
]}
]
, {t, 0, 10, 
        10/\((frames - 1)\)}]\)

Problem 4: Springs Mounted on Top of a Carriage (corrected; 20 Sept. 2006)

Solution

Clear["Global`*"]

Part a

T=(1/2) m (((xi'[t])^2)+((yi'[t])^2));

V = 
        (1/2) K (((-r0 + xb[t]))^2) + (1/2) k yb[t]^2;

change = {
xi -> (((xb[#]) 
                      Cos[ω #] - yb[#] Sin[ω #]) &), 
          yi -> (((xb[#]) Sin[ω #] + 
                    yb[#] Cos[ω #]) &)};

L =T-V  /.change  //FullSimplify

eqmotion  =  
      EulerEquations[L, {xb[t], yb[t]}, t] //  
        Solve[#,  {xb''[t], yb''[t]}] &

Part b (corrected; 20 Sept. 2006)

values = {K ->  1, k -> 1 , m -> 1, 
          r0 -> 1};

{eq1, 
        eq2} = (eqmotion[[1]] /. values /. 
          Rules -> Equal)

eq3 = xb''[t] + 
          I yb''[t]  == (xb''[t] + I yb''[t] /. 
              eq1 /. eq2)

eq4 = 
      eq3 /. xb -> ((u[#] - I  yb[#]) &) // 
        Simplify

uSol = u[
                t] /. DSolve[{eq4,  
                    u[0] == 1, u'[0] == 0 } ,  u[t], 
                  t][[1]] //   ExpToTrig // Simplify // 
        Expand

xSol =   uSol /.  
          Complex[a_, b_] -> 
            Complex[a, 0]//Simplify

ySol = (uSol - xSol  // Together) /. 
            I -> 1 // Simplify

ySol = (uSol - xSol  // Together)/I  // 
          Simplify

wvalues = {0.2, 0.4, 0.6, 0.8};

\(\(listb\  = \(ParametricPlot[\ {xSol, ySol} /. ω -> #\  // 
                Evaluate, {t, 0, 100}, 
              PlotPoints -> 100, 
AspectRatio -> 
                1, 
AxesLabel -> {"xb", "yb"}, 
Ticks -> 
                False, 
DisplayFunction -> 
                Identity] &\)\ \  /@ \ wvalues;\)\)

Show[GraphicsArray[Partition[listb , 2]]];

change

eq2 = {xi[t], yi[t]} /. change /. xb[t] -> xSol /. 
          yb[t] -> ySol  // Simplify

\(\(listi\  = \(ParametricPlot[
              eq2 /. ω -> # // Evaluate, {t, 0, 
                100}, 
AxesLabel -> {"xi", "yi"}, 
AspectRatio -> 
                1, 
Ticks -> 
                False, 
DisplayFunction -> 
                Identity] &\) /@ \ wvalues\ ;\)\)

Show[GraphicsArray[Partition[listi , 2]]];

Part c

listpb  = ParametricPlot[ {ySol, 
                    D[ySol, t]}   /. ω -> # // Evaluate, {t, 0, 
                100}, 
AspectRatio -> 
                1, 
AxesLabel -> {"yb", " \
yb'"}, 
Ticks -> 
                False, 
DisplayFunction -> 
                Identity] &  /@   wvalues;

Show[GraphicsArray[Partition[listpb , 2]]];

listpi  = ParametricPlot[ {eq2[[2]], 
                    D[eq2[[2]], t]}   /. ω -> # // 
                Evaluate, {t, 0, 
                100},  
AspectRatio -> 
                1, 
AxesLabel -> {"yb", " \
yb'"}, 
Ticks -> 
                False, 
DisplayFunction -> 
                Identity] &  /@   wvalues;

Show[GraphicsArray[Partition[listpi , 2]]];

Problem 5: Mass Falling Through a Hole in a Table

Solution

Clear["Global`*"]

Part a

T1=(1/2) m1 (((r'[t])^2)+r[t]^2 ((ϕ'[t])^2));
V1 = 0;

zRule = {z -> ((len - r[#]) &)};

T2=(1/2) m2 ((z'[t])^2)/.zRule;
V2 = -m2 g z[t] /. zRule;

L=(T1+T2)-(V1+V2)

eq1= EulerEquations[L, {r[t], ϕ[t]}, 
        t]

eq2 =  
        Solve[eq1, {ϕ''[t], r''[t]}] // 
          Flatten;
eq2 // ColumnForm

FirstIntegrals[L, {r[t], ϕ[t]}, t] // 
      ColumnForm

eqϕ= {ϕ'[t]->(h/r)[t]^2,ϕ''[t]->-((2 h r'[t])/(r[t]^3))};

eqR =  eq2[[1]] /. eqϕ // Simplify // 
        ExpandAll

eqMotion = {eqϕ[[1]], eqR} /. Rule -> Equal // 
          Flatten;
eqMotion // ColumnForm

Part b

energy = T1 + T2 + V1 + V2 /. eqϕ // Simplify // 
        ExpandAll

Effpot  = (+g  len  m2 + energy) /. {r'[t] -> 0, 
            r[t] -> r} // Apart

min = 
      D[Effpot, r ] == 0 // Solve[#, r ][[1]] &

values = {h -> 1, m1 -> 1, m2 -> 1, 
          g -> 1 };

Plot[Effpot /. values // Evaluate, {r , 0.2, 3}];

ContourPlot[ energy /. values /. len -> 10  // 
          Evaluate, 
{r[t], 0.2, 6}, {r'[t], -1, 
          1}, 
ContourShading -> 
          True, 
ColorFunction -> 
          Hue  ];

Part c

values = {h -> 1, m1 -> 1, m2 -> 1, 
          g -> 1};

ndsol  = NDSolve[{ eqMotion /. values , ϕ[0] == 0, 
                r[0] ==  2, r'[0] == 0} // 
              Flatten, {ϕ[t], r[t]}, {t, 0, 
              70}][[1]] ;

ϕ1[t_]= ϕ[t] /.ndsol; 
r1[t_]= r[t] /.ndsol;

ParametricPlot[ 
        r1[t] {Cos[ϕ1[t]], Sin[ ϕ1[t]]}   // Evaluate, {t, 0, 
          70}, 
AspectRatio -> 
          Automatic, 
Epilog -> {AbsoluteThickness[3], 
            Circle[{0, 0}, 2], 
            Circle[{0, 0},  .55]}
];

Problem 6: Spring Pendulum

Solution

Clear["Global`*"]

Part a

lsrule = 
        lst -> ((Sqrt[x[#1]^2 + y[#1]^2]) &);

T= (1/2) m (((x'[t])^2)+((y'[t])^2));
Vsp =  
        (1/2) k (((-L0 + lst[t]))^2);
Vgr = - m g x[t]

polarrule = {x -> ((r[#] Cos[θ[#]]) &), 
          y -> ((r[#] Sin[θ[#]]) &)};

L = T - Vsp - Vgr /. lsrule /. polarrule // 
          Simplify//PowerExpand

eqMotion = ((  EulerEquations[L, {r[t], θ[t]}, 
                    t]
 // ( Solve[#, {r''[
                              t], θ''[
                              t]}][[1]] & )
) /. 
              Rule -> Equal 
 // 
            ExpandAll
);
eqMotion  // ColumnForm

Part b

vertical =   {θ -> (0 &)};

vertMotion = eqMotion /. 
          Rule -> Equal/.vertical

vertSol = DSolve[vertMotion[[1]] ,  r[t], t] // Flatten // 
        ExpandAll

vertSol /. C[1] -> 0 /. C[2] -> 0 // Flatten // 
      ExpandAll

Part c

valuesc = {L0 -> 40/9, n -> 1, k -> 4,  g -> 1, 
          m -> 1};

initialc = {r[0] == (L0 +  g m/k) + 1 , r'[
                0] ==  .1, θ[0] == 
              1/2   π, θ'[0] ==   .0} /. 
          valuesc;

neqc = {eqMotion, initialc} /. Rule -> Equal  /. valuesc // 
          Flatten ;

nSolc = 
        NDSolve[neqc, {r[t], θ[t]}, {t, 0, 110}, 
            MaxSteps -> 2000 ] // Flatten;

rc[t_]= r[t]/.nSolc; 
θc[t_]= θ[t]/.nSolc;

Plot[ {rc[t], θc[t]  } // Evaluate, {t, 0, 
          110}, 
PlotStyle -> {Hue[ .7], 
            Hue[ .9]}  , 
Ticks -> {Automatic,  π Range[-1, 2, 1  ] } , 
Epilog -> 
          Map[Line,  {{0, #  π}, {125, #  π} } &  /@   Range[-1, 2]] 
];

ParametricPlot[ 
        rc[t] {Sin[θc[t]], -Cos[θc[t]]}  // Evaluate, {t, 
          0, 110}, 
PlotStyle -> 
          Hue[0.9],  
AspectRatio -> 
          1, 
PlotRange -> {Automatic, {-8, 
              8}},  
Epilog -> { Circle[{0, 
                    0}, #] & /@ {6, 4, 
                2} }
];

pt1 = ParametricPlot[ {#[t], #'[t]  } // Evaluate, {t, 0, 
                110}, 
PlotStyle -> {Hue[ .7], 
                  Hue[ .9]}, 
Ticks -> 
                False, 
DisplayFunction -> 
                Identity 
] &  /@   {rc, θc};

    
Show[GraphicsArray[pt1]];

Part d

valuesd = {L0 -> 5, n -> 1, k -> 2,  g -> 1, 
          m -> 1};

initiald = {r[0] == (L0 +  g m/k +  .01)  , r'[
                0] == 
              0, θ[0] == π  -  .07, θ'[
                0] ==  0} /. valuesd ;

neqd = {eqMotion, initiald} /. Rule -> Equal  /. valuesd // 
          Flatten;

nSold = 
        NDSolve[neqd, {r[t], θ[t]}, {t, 0, 110}, 
            MaxSteps -> 2000] // Flatten;

rd[t_]= r[t] /. nSold;
θd[t_]= θ[t] /. nSold;

Plot[ {rd[t], θd[t]  } // Evaluate, {t, 0, 
          110}, 
PlotStyle -> {Hue[ .7], 
            Hue[ .9]} , 
Ticks -> {Automatic,  π Range[-5, 2, 1  ] } , 
Epilog -> 
          Map[Line, {{0, #  π}, {125, #  π} } &   /@   Range[-5, 2]] 
];

ParametricPlot[ 
        rd[t] {Sin[θd[t]], -Cos[θd[t]]} // Evaluate, {t, 0, 
          110}, 
PlotStyle -> 
          Hue[ .9],  
PlotRange -> 
          All, 
AspectRatio -> 
          Automatic,   
Epilog -> { Circle[{0, 
                    0}, #] & /@ {8, 6, 
                4}}
];

pt1 = 
ParametricPlot[ {#[t], #'[t]  } // 
                Evaluate, {t, 0, 
                110}, 
PlotStyle -> {Hue[ .7], 
                  Hue[ .9]} , 
Ticks -> 
                False, 
DisplayFunction -> 
                Identity 
] &  /@  {rd, θd};

    
Show[GraphicsArray[pt1]];

Part e

Contains correction on:  
    2nd Edition:Chpt.6,Problem 6.2.6 p.334:Spring Pendulm. 
(FIO, 17 Jan 2003)  

All values should be consistently either {rc,θc} or {rd,θd}

line[t_]:=
Show[ 
Graphics[
{{Line[{{0, 0}, 
                      rc[t] {Sin[θc[t]], -Cos[θc[t]]}}],
{PointSize[0.1],Point[ rc[t]{Sin[θc[
                                t]], -Cos[θc[t]]}]}}},
PlotRange -> {{-10, 
                  10}, {-13, 13}}
]
]

frames = 2;
movie = line[#] &  /@  
          Range[0, 30, 30/( frames - 1)];

6.3 Lagrangian Problems with Lagrange Multipliers

Overview of Nonholonomic Constraints and Lagrange Multipliers

Lag[qi_,L_ ,Qi_,t_:t] := 
D[ D[L,(D[qi[t],t])],t]-D[L,qi[t]]==Qi;
  
Lag[qi_List,L_,Qi_List,t_:t] :=
(D[ D[L,(D[qi[[#]][t],t])],t]-D[L,qi[[#]][t]]==Qi[[#]]
)&  /@  Range[1,Length[qi]];

Protect[Lag];

Problem 1: Atwood Machine

Remarks and outline:

Solution

Clear["Global`*"]

T=(1/2) Sum[m[i] (((x[i])'[t])^2),{i, 1, 2}];
V = (-g) Sum[
            m[i] x[i][t], {i, 1, 2}];

L= T-V

f[1]=x[1][t]+x[2][t]-Len;

Qeq = D[f[1], x[#][t]] λ[1] &  /@  {1, 
          2}

eq2=Lag[{x[1], x[2]}, L, Qeq]

f[1] == 0

eq3 = D[f[1] == 0, {t, 2}]

eq4 = ((Solve[{eq2,  eq3} // 
                      Flatten , { λ[1], (x[2])''[
                        t], (x[1])''[t] }, t]  // 
                  Simplify
 // 
                Flatten 
) /. {Rule -> 
                Equal}
);
eq4 // ColumnForm

initial = {x[1][0] == x0, (x[1])'[0] == 
            v0};

sol1 =  DSolve[{eq4[[{3}]], initial}  // 
              Flatten, {(x[1])[t]}, t] // 
          FullSimplify //Flatten

sol2 = {x[2][t] -> Len - x[1][t] /. sol1}

Problem 2: Hoop Rolling on an Incline

Remarks and outline

Solution

Clear["Global`*"]

Part a

T=(1/2) m (((x'[t])^2)+(r^2) ((θ'[t])^2));
V = 
        m g (Len - x[t]) Sin[ϕ];
L = 
      T - V

eq1 = r θ[t] - x[t] == 0;

constraintForce = 
{Q[
            1] -> λ D[eq1[[1]], x[t]], 
Q[
            2] -> λ D[eq1[[1]], θ[t]]}

{eqX, eqθ} = 
        Lag[{x, θ}, L, {Q[1], Q[2]}] //. 
          constraintForce;
{eqX, eqθ} // ColumnForm

Part b

initial = {x[0] == x0, x'[0] == 
            v0, θ[0] == x0/r, θ'[0] == 
            v0/r};

xtSol = DSolve[{eqX, eqθ, initial} // 
            Flatten, {x[t], θ[t]}, t][[1]]

λSol = Solve[eq1 /. xtSol, λ]

xtSol2= xtSol /.λSol //ExpandAll //Flatten

Part c

values = {x0 -> 0, v0 -> 0, 
          g -> -0.1, ϕ -> π/4, 
          r -> 2.0};

posx = ( x[t] Sin[ϕ]+ r Sin[θ[t]]);
posy = ( x[t] Cos[ϕ]+ r Cos[θ[t]]);

point={posx,posy} /.xtSol2 //.values;

ParametricPlot[point // Evaluate, {t, 0, 70}, 
        AspectRatio -> 1 ];

Problem 3: Sphere Rolling on a Fixed Sphere

Contains correction on:  
    2nd Edition: Chpt.6, Problem 6.3.3 p.174: Sphere Rolling on a Fixed Sphere.
(FIO, 17 Jan 2003)  

The figure in the first edition was incorrect. This is fixed in the 2nd edition. The second angle (theta-2) is measured from the contact point of the two spheres. Although this definition may be less conventional, it simplifies our particular problem.

Remarks and outline

Solution

Clear["Global`*"]

Part a

Tcm=(1/2) m (((r'[t])^2)+r[t]^2 (((θ1)'[t])^2));

Trot=(1/2) I1 ((((θ1)'[t]+(θ2)'[t]))^2);

T=Tcm+Trot;

V= m g r[t] Cos[θ1[t]];

L=T-V

constraints= {b θ1[t] - a θ2[t]==0, r[t]-a-b==0 }

constraintRule= DSolve[constraints , {θ2, r}, 
        t]

vars = {r[t], θ1[t], θ2[t]};

Qeq = Sum[λ[j]*D[constraints[[j, 1]], vars[[#]]], {j, 
              1, 2}] &  /@  {1, 2, 3}

eqMotion =  
        Lag[{r, θ1, θ2}, L, Qeq];
eqMotion // ColumnForm

Part b

{λ2Sol, θSol, λ1Sol} = (Solve[
                eqMotion //. constraintRule // 
                  Flatten, {θ1''[t], λ[1], λ[
                    2]}]
 /. 
              Rule -> Equal
 // 
            Flatten
);

    
{λ2Sol, θSol, λ1Sol} // 
      ColumnForm

eqθ=( θSol[[1]] - θSol[[2]])

intθ = Integrate[
              eqθ  θ1'[t], {t, 0, 
                t}]  //. {θ1[0] -> 0, θ1'[0] -> 
                0} // Simplify // Apart

solθ = Solve[
            intθ == 
              0, {θ1'[
                t]}][[1]]  /. (Sqrt[-1 + Cos[x_]]) -> I (Sqrt[1 - Cos[x]])

Part c

eq1 = λ2Sol /. λ[2] -> 0 //. 
        solθ

sol2=Solve[eq1,Cos[θ1[t]]] //Flatten

sol3 = Solve[sol2 //. {Rule -> Equal}, θ1[t], 
          InverseFunctions -> 
            True]//Flatten

(θ1[t]/Degree) //. 
          sol3[[2]] //. {I1 -> ((2 (a^2) m)/5)} // N

6.4 Hamiltonian Problems

Overview of Hamilton's Equations

Hamilton[L_, xList_List, pList_List, 
        t_:  t] := 
Module[{xx, vv, pp, sol, ham, eqp, eqx, 
          eqs}, 
xx = 
          Through[xList[t]]; 
vv = 
          D[xx, t]; 
pp = 
          Through[pList[t]]; 
sol = 
          Solve[(D[L, #] & /@ vv) == pp, vv] // 
            Flatten; 
ham = pp . vv - L //. sol // 
              Simplify // Expand; 
eqp = 
          D[pp, t] == -Map[D[ham, #] &, xx] // 
            Thread; 
eqx = 
          D[xx, t] == +Map[D[ham, #] &, pp] // 
            Thread; 
eqs = 
          Join[eqp, eqx]; 
Return[{sol, ham, 
            eqs}]
]

Clear[L, x, V]

L=-V[x[t]]+(1/2) m ((x'[t])^2);

{pRule, hamiltonian, eqMotion} = 
        Hamilton[L, {x}, {px}];
pRule
hamiltonian
eqMotion

L=-V[x[t], y[t]]+(1/2) m (((x'[t])^2)+((y'[t])^2));
Hamilton[
      L, {x, y}, {px, py}]

Protect[Hamilton];

Problem 1: Harmonic Oscillator

Remarks and outline

Solution

Clear["Global`*"]

T=(1/2) m ((x'[t])^2);
 V = (1/2) k x[t]^2;

L = T - V

eq1 = p[t] == D[L, x'[t]]

xpRule = Solve[eq1, x'[t]][[1]]

H = (p[t] x'[t] - L) //. xpRule

eq2 = 
{p'[
            t] == -D[H, x[t]], 
x'[
            t] == +D[H, p[t]]}

Hamilton[L, {x}, {p}] // ColumnForm

initial = {p[0] == p0, x[0] == x0};

dsol = DSolve[eq2~Join~initial, {p[t], x[t]}, t] // Flatten  // 
        FullSimplify

Problem 2: Nonlinear Oscillator

Clear["Global`*"]

Remarks and Outline

firstDiffSeries[eqin_, {z_, p_}, initRule_List, order_, 
        t_:  t] := 
Module[{eq, zSer, pSer, eSer, vars, 
          sol, zpsol}, 
eq = 
          Table[eqin[[i, 1]] - eqin[[i, 2]] == 0, {i, 
              2}]; 
zSer = 
          Series[z[t], {t, 0, order}]; 
pSer = 
          Series[p[t], {t, 0, order}]; 
eSer = 
          eq //. ({z[t] -> zSer, p[t] -> pSer}) // 
            ExpandAll; 
vars = (Table[{D[z[t], {t, j}], 
                    D[p[t], {t, j}]}, {j, 1, order + 1}] /. t -> 0) // 
            Flatten;      sol = 
          Solve[eSer, vars] // 
            Flatten; 
zpSol = {zSer, pSer} //. sol //. 
            initRule; 
Return[
          zpSol]
]

eq1={(w^2) x[t]+p'[t]==0,p[t]==x'[t]};
firstDiffSeries[
        eq1, {x, p}, {x[0] -> x0, p[0] -> p0}, 4, t] // 
      ColumnForm

firstOrderPert[eq1_, x_, {x0_, v0_}, eps_:  eps, 
        t_:  t] := 
Module[{xRule, eq2, eqEps, initCon0, 
          initCon1, dsol1, 
          dsol2}, 
xRule = {x -> (x[0][#] + 
                    eps x[1][#] &)}; 
eq2 = 
          eq1[[1]] - eq1[[2]] /. xRule; 
eqEps = 
          CoefficientList[eq2, eps] == 0 // 
            Thread; 
initCon0 = {x[0][0] == 
              x0, (x[0])'[0] == 
              v0}; 
initCon1 = {x[1][0] == 
              0, (x[1])'[0] == 
              0}; 
dsol1 = (DSolve[
                Join[{eqEps[[1]]}, initCon0], {x[0][t]}, t] // 
              Flatten); 
dsol2 = (DSolve[
                Join[{eqEps[[2]]}, initCon1] //. dsol1, {x[1][t]}, 
                t] // Flatten); 
Return[x[t] /. 
              xRule //. Join[dsol1, dsol2]]
]

eq1=x''[t]==-(k^2) x[t] + eps x[t]^2;

firstOrderPert[eq1, x, {x0, 0}, eps, t] // Simplify  // 
      Collect[#, eps] &

Solution

Part a

T=(1/2) m ((x'[t])^2);
V = 
        (1/2) k x[t]^2 - (1/4) ϵ x[t]^4;

L  = T - V;

{xprule, hamiltonian, eqMotion} = 
        Hamilton[L, {x}, {p}];
{xprule, hamiltonian, eqMotion} // ColumnForm

Part b

values = {k ->  1,  ϵ -> 1/2, 
          m -> 1};

pt1 = 
        Plot[V  /. values // Evaluate, {x[t], -2 , 2 }, 
          DisplayFunction -> Identity];

pt2 =  
        ContourPlot[
          hamiltonian  /. values // Evaluate, {x[t], -2, 
            2 }, 
{p[t], -1, 
            1}, 
ColorFunction -> 
            Hue, 
Contours -> 
            15, 
DisplayFunction -> 
            Identity] ;

Show[GraphicsArray[{pt2, pt1}]];

Part c

solt = firstDiffSeries[
            eqMotion, {x, p}, {x[0] -> x0, p[0] -> 0}, 6] // 
          Normal  // Simplify

eq1 = (eqMotion // First) /. {p'[t] -> 
            m x''[t]}

solϵ =firstOrderPert[eq1, x, {x0, 0}, ϵ]//Collect[#, {x0, ϵ,   Cos[x_] }] & //Map[FullSimplify, #, 1] &

test1 = solt[[1]] /. ϵ^n_ /; n > 1 -> 0 // Apart // 
        Collect[#, t] &

test2 =  
      Series[solϵ, {t, 0, 6}] // Normal

test1 == test2

values = {x0 -> 1.2 , p0 -> 0, k -> 1, 
          m -> 1, ϵ -> 1/2};

nsol = 
        NDSolve[{eqMotion /. values, p[0] == 0, 
              x[0] == 1.2  }, {p[t], x[t]}, {t, 0 , 30 }] // 
          Flatten;

Plot[{ solϵ , solt[[1]], x[t]} /. nsol /. 
            values    // Evaluate, {t, 0, 
          30}, 
PlotStyle -> {Dashing[{0.01}], 
            Dashing[{0.05}],  
            Thickness[0.005]}, 
GridLines -> 
          Automatic] ;

Problem 3: Cylindrical Coordinates

Clear["Global`*"]

Remarks and outline

firstDiffSeries[eqin_, {z_, p_}, initRule_List, order_, 
        t_:  t] := 
Module[{eq, zSer, pSer, eSer, vars, 
          sol, zpsol}, 
eq = 
          Table[eqin[[i, 1]] - eqin[[i, 2]] == 0, {i, 
              2}]; 
zSer = 
          Series[z[t], {t, 0, order}]; 
pSer = 
          Series[p[t], {t, 0, order}]; 
eSer = 
          eq //. ({z[t] -> zSer, p[t] -> pSer}) // 
            ExpandAll; 
vars = (Table[{D[z[t], {t, j}], 
                    D[p[t], {t, j}]}, {j, 1, order + 1}] /. t -> 0) // 
            Flatten;      sol = 
          Solve[eSer, vars] // 
            Flatten; 
zpSol = {zSer, pSer} //. sol //. 
            initRule; 
Return[
          zpSol]
]

Required packages

Needs["Calculus`VectorAnalysis`"]

Solution

Part a

rRule = 
      rVector -> (Evaluate[
              CoordinatesToCartesian[{r[#], θ[#], z[#]}, 
                Cylindrical]] &)

T = ((1/2) m rVector'[t] . rVector'[t]) //. rRule // 
        Simplify
L = T - V[r[t], θ[t], z[t]] // Simplify

{pRule, hamiltonian, eqMotion} = 
        Hamilton[L, {r, θ, z}, {pr, pθ, pz}];

pRule

hamiltonian

eqMotion // ColumnForm

Part b

vRule = {V -> (-(k/(Sqrt[#(1^2) + #(3^2)])) &)};
V[r[t], θ[t], z[t]] //. vRule

zRule = {z -> (0 &), pz -> (0 &)};

eq1 = eqMotion //. vRule /. zRule // 
          PowerExpand;
eq1 // ColumnForm

eq2 = eq1 //. pθ -> (m h &)

rEq = eq2[[{1, 4}]];
rEq // ColumnForm

Part c

initial = {r[0] -> r0, pr[0] -> 0};

pert =  firstDiffSeries[rEq, {r, pr}, initial, 4, 
          t]//Simplify

Problem 4: Swinging Atwood Machine

Solution

Clear["Global`*"]

Part a

T = (m/2) {x1'[t], y1'[t]} .  {x1'[t], y1'[
                  t]} +  (M/2) x2'[t]^2;
 V = -g  m x1[t] - g   M x2[t];

change = 
{x1 -> (r[#] Cos[θ[#]] &), 
y1 -> (r[#] Sin[θ[#]] &)};

L0Rule = x2 -> (L0 - r[#] &);

L=T-V /.L0Rule /.change //Simplify

{xpRule, ham, eqMotion} = 
        Hamilton[L , {θ, r}, {pθ, pr}] // 
          Simplify;

xpRule

ham

eqMotion // ColumnForm

eqMotion /. {pθ[t] -> 0,  θ[t] -> 
            0} /. r ->(x1[#] & )

Part b

values = {m -> 1, M -> 1.05, g -> 1};

initial = {r[0] == 10, 
          pr[0] == 0, θ[0] == π/12, 
          pθ[0] == 0};

eq1 = {eqMotion, initial} /. values // Flatten;

nSol = 
        NDSolve[eq1, {r , pr, θ , pθ}, {t, 0, 250 }  ] // 
          Flatten ;

{r1[t_], θ1[t_], pr1[t_], 
          pθ1[t_]} = {r[t], θ[t], pr[t],  pθ[t]} /. 
          nSol;

Plot[{ r1[t],  θ1[t]} // Evaluate, {t, 0, 
          250}, 
PlotStyle -> {Dashing[{0.01}], 
            Thickness[0.008]} ];

ParametricPlot[{ 
            r1[t] Sin[θ1[t]], -r1[t] Cos[θ1[t]]} // 
          Evaluate, {t, 0, 80}];

pt1 = 
ParametricPlot[# // Evaluate, {t, 0, 
                250}, 
Ticks -> 
                False, 
DisplayFunction -> 
                Identity
] &  /@   {{ r1[t]  ,  
              pr1[t]  }, { θ1[t]  ,  
              pθ1[t]  }};

Show[GraphicsArray[pt1]];

Problem 5: Spherical Pendulum

Required Packages

Needs["Graphics`Shapes`"]

Solution

Clear["Global`*"]

Part a

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

change = 
{x  -> (L0 Cos[ϕ[#]]  Sin[ θ[#]] &), 
y  -> (L0 Sin[ϕ[#]]  Sin[ θ[#]] &), 
z  -> (L0 Cos[θ[#]] &)};

ψRule =  θ -> (π - ψ[#] &);

L = T - V /. change /. ψRule // Simplify

eq1 =  {prule, ham, eqMotion} = 
          Hamilton[L, {ψ, ϕ}, {pψ, pϕ}, 
            t];

prule

ham

eqMotion  // ColumnForm

eqMotion // ColumnForm

pϕRule = {pϕ -> (m h &)};

eqMotion /. pϕRule // ColumnForm

eqMotion  /.  pϕRule /. h -> 0 // 
        ExpandAll//Flatten

Solve[{eqMotion[[1]],  
                  D[eqMotion[[3]], t]} , ψ''[t]]  /. 
              Rule -> Equal /.  pϕRule /. h -> 0 // 
        ExpandAll//Flatten

Part b

eqMotion

values = {L0 -> 1, g -> 1, m -> 1};

initial = {ϕ[0] == 0, 
          pϕ[0] == 0.1, ψ[0] == π/2, 
          pψ[0] == 0};

nSol  =  
        NDSolve[  {eqMotion, initial} /. values // 
              Flatten , {ψ, ϕ, pψ, pϕ},  {t,  0,  
              25}]  // Flatten ;

{ψ1[t_] , ϕ1[t_] } =   { ψ[t],   ϕ[t] } /. 
          nSol;

pt1 = 
ParametricPlot3D[
{x[t], y[t], 
                        z[t], Thickness[ .01]}  /. change /. ψRule /. 
                  values /. ψ -> ψ1 /. ϕ -> ϕ1  // Evaluate, 
{t, 0, 
            25}, 
DisplayFunction -> 
            Identity
];

Needs["Graphics`Shapes`"]

sphere =  Graphics3D[Sphere[0.96]]  ;

pt1 = 
        Show[ sphere, pt1, Boxed -> False , Shading -> False, 
          ViewPoint -> {1, 1, 0}];

Show[pt1 /. Polygon[x_] -> {Hue[0.8], Line[x]},  
        ViewPoint -> {1, 0, 1/2}];

6.5 Hamilton-Jacobi Problems

Overview

Clear["Global`*"]

HamiltonJacobi[L_, xList_List, pList__List, wList_List, a_, 
        t_:  t] := 
Module[{xx, vv, pp, sol, ham, ww, 
          hamPrincFunc, pNewRule, HJeq}, 
xx = 
          Map[#[t] &, xList]; 
vv = 
          Map[#'[t] &, xList]; 
pp = 
          Map[#[t] &, pList]; 
sol = 
          Solve[(D[L, #] & /@ vv) == pp, vv] // 
            Flatten; 
ham = (pp . vv - L) //. sol // 
              Simplify // Expand; 
ww = 
          Map[#[t] &, 
            wList, {2}]; 
hamPrincFunc = (Plus @@ ww) - 
            a t; 
pNewRule = {pp -> (D[
                          hamPrincFunc, #] & /@ xx) // Thread} // 
            Flatten; 
HJeq = (ham //. pNewRule) - 
              a == 0;        
Return[{sol, 
            hamPrincFunc, HJeq}]
]

L=+(1/2) m ((x'[t])^2)-V[x[t]];
{pRule, hamPrincFunc, eqHamJac} = 
        HamiltonJacobi[L, {x}, {p}, {W[x]}, a];

pRule

hamPrincFunc

eqHamJac  // ColumnForm

L=+(1/2) m (((x'[t])^2)+((y'[t])^2))-V[x[t], y[t]];
HamiltonJacobi[
      L, {x, y}, {px, py}, {W1[x], W2[y]}, a]

Protect[ HamiltonJacobi ];

Problem 1: Harmonic Oscillator

Solution

Clear["Global`*"]

Part a

T=(1/2) m ((x'[t])^2);
V = 
        (1/2) k x[t]^2;
L = T - V

{xpRule, ham, eqMotion} = 
        Hamilton[L, {x}, {p}] //. {x_[t] -> x};

xpRule

ham

eqMotion

pRule = {p -> D[S[x, t], x]};

eqHJ = (ham + D[S[x,t],t]) ==0  //. pRule

sRule = {S -> (W[#1] - P #2 &)};

eqW= eqHJ //.sRule

Part b

dWRule= Solve[eqW,W'[x]][[2]]

wRule = 
      W[x] ->   (Sqrt[2  m P])  Integrate[(Sqrt[1 - (x^2) λ]), 
              x] /. λ -> (k/(2 P))

eq1= Q== D[S[x,t] /.sRule /.wRule,P] //Simplify

xSol= Solve[eq1,x ] //Simplify //PowerExpand //Simplify //Flatten

eq2 = xpRule /. Rule -> Equal // Solve[#, p] &  // 
        Flatten

eq2 /. x' -> D[x /. xSol, t]

Problem 2: Particle in a Constant Gravity Field

Solution

Clear["Global`*"]

Part a

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

{prule, S, Chareq} = 
      HamiltonJacobi[L, {x, z}, {px, pz}, {Wx[x] , Wz[z]}, P] // 
        ExpandAll

eq1a = 
      Chareq[[1, {1, 2, 4}]] == A  /. 
        z[t] -> z
eq1b = 
      Chareq[[1, 3]] == -A    /. x[t] -> x

Part b

Wxsol = DSolve[eq1a, Wz[z], z][[2]]
Wzsol = DSolve[eq1b, Wx[x], x][[2]]

Ssol = -P t + Wx[x] + Wz[z] /.  
          Wxsol/.Wzsol

Qeq = {Qx, 
              Qz} == ( D[Ssol, #] & /@ {A, P} // Flatten) // 
          Thread //Simplify

Sol = Solve[Qeq, {x, z}] /. {x -> x[t], z -> z[t]} /. 
          Rule -> Equal // Flatten

\(\(\(\(Solve[{Sol, \((Sol /. t -> 0)\), \((D[Sol, t] /. 
                  t -> 0)\)} // Flatten, {x[t], z[t]}, {A, P, Qx, Qz}] // 
        Flatten\) // ExpandAll\)\(
\)
    \)\)

Problem 3:  Kepler's Problem and Hamilton-Jacobi Equations

Remarks and outline

Required package

Needs["Calculus`VectorAnalysis`"]

Solution

Clear["Global`*"]

Part a

Needs["Calculus`VectorAnalysis`"]

rRule = 
      rvector -> (Evaluate[
              CoordinatesToCartesian[{r[#], θ[#], z[#]}, 
                Cylindrical]] &)

T = (m/2) rvector'[t] . rvector'[t] //. 
              rRule  //. {z'[t] -> 0}  // 
          Simplify;
V = -(k/r)[t];

L = T - V

{momRule, hamPF, eqMotion} = 
        HamiltonJacobi[
          L, {r, θ}, {pr, pθ}, {Wr[r], Wθ[θ]}, 
          En];

    
{momRule, hamPF, eqMotion} // ColumnForm

wθrule = DSolve[{Wθ'[θ] == +el, 
            Wθ[0] == 
              0}, {Wθ}, θ][[1]]

eq1 = (eqMotion //. 
                wθrule // Solve[#, Wr'[
                      r[t]]][[1]] &) /.   r[t] -> r  // 
        PowerExpand

eq2 = D[eq1[[1, 2]], el]

eq3 = 
      Integrate[ (1/(r (Sqrt[-1 + a r + b (r^2)]))), 
          r ] /.  {a -> ((2 k m)/(el^2)), 
          b -> ((2 En m)/(el^2))}

Qeq =  
      Q == θ +eq3

Part b

ueq = Qeq /. r -> (1/u)

usol = u  /. Solve[ueq, u][[2]]  /. {Q -> π/
                    2 } // Simplify // PowerExpand // 
        Simplify

eq4 = 
      usol == ((1 + e Cos[θ])/(a (1 - (e^2))))

eq5 =  eq4 /. {θ -> {0, π/2 }}  // Thread

Solve[eq5 ,{a,e}]

6.6 Exercises

Exercise 1: Atwood Machine

Exercise 2:  Double Atwood machine

Exercise 3: Spherical Pendulum

Exercise 4: Double Pendulum

Exercise 5:  Spring Pendulum

Exercise 6:  Hamilton-Jacobi in parabolic coordinates

Exercise 7:  How Hamilton Works

Exercise 8:  How HamiltonJacobi Works

Exercise 9:  How firstDiffSeries Works

Exercise 10: How firstOrderPert Works