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.1 Introduction

4.1.3 Mathematica Commands

Example L1:  Lagrange's equation for a particle in a potential V[x]

L1 = m/2 x'[t]^2  -  V[x[t]]

leq1  = { Lag[x, L1]}

L2 = m/2 (x'[t]^2  +  y'[t]^2)  -  V[x[t], y[t]]

leq2  = Lag[{x, y}, L2]

Example H1:  Hamilton's equations of motion in one dimension.

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

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

pRule

ham

eqMotion

Example H2:  Hamilton's equations of motion in two dimensions.

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

{pRule, ham, eqMotion} = Hamilton[L, {x, y}, {px, py}];

pRule

ham

eqMotion

How Hamilton works.

Version of Hamilton with no local variables:

\(\(xHamilton[L_, xList_List, pList_List, t_:  t] := \ 
\ \ 
      Module[{}, 
\ \ \ \ 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; 

        \ \ \ \ 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}]

        \ \ \ ]
\)\)

1-D Lagrangian:

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

Hamilton results:

xHamilton[L, {x}, {px}]

Representation of position, velocity, and momentum:

{xx, vv, pp}

Equation relating velocity and momentum:

(D[L, #]& /@ vv)  ==  pp

Solution of this converts velocity to momentum:

sol

Hamiltonian defined by H = mom × vel  -  L :

ham

Equations of motion given by dp/dt = -dH/dx and dx/dt = +dH/dp:

{eqp, eqx}

Example HJ1:  One-dimensional particle in a potential V[x].

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

\(\({pRule, hamPrincFunc, eqHamJac} = 
	\ \ 
      HamiltonJacobi[L, {x}, {p}, {W[x]}, a]; \)\)

pRule

hamPrincFunc

eqHamJac

Example HJ2:  Two-dimensional particle in a potential V[x].

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

\(\({pRule, hamPrincFunc, eqHamJac} = 
	\ \ 
      HamiltonJacobi[L, {x, y}, {px, py}, {W1[x], W2[y]}, a]; \)\)

pRule

hamPrincFunc

eqHamJac

How HamiltonJacobi works.

Version of HamiltonJacobi with no local variables:

\(xHamiltonJacobi[L_, xList_List, pList__List, wList_List, a_, t_:  t] := 
      \ 
\ \ Module[{}, 
\ \ \ \ 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}]
\ \ \ ]\)

1-D Lagrangian:

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

Hamilton results:

xHamiltonJacobi[L, {x}, {p}, {W[x]}, a]

Representation of position, velocity, and momentum are the same as in Hamilton.  The W list is:

ww

The Hamilton's Principle function is:

hamPrincFunc

Rule to eliminate momentum in terms of W'[x[t]]

pNewRule

The Hamilton-Jacobi equation given by H - a with this momentum substitution:

HJeq

Example DS1:  Expansion of harmonic oscillator using a series expansion of two first-order equations.

eq1  =  {p'[t]  + w^2 x[t]  ==  0,  p[t] == x'[t]}

{xser, pser} = 
      firstDiffSeries[eq1, {x, p}, {x[0] -> x0, p[0] -> p0}, 4, t];

xser

pser

How firstDiffSeries works.

Version of firstDiffSeries with no local variables:

\(xfirstDiffSeries[eqin_, {z_, p_}, initRule_List, order_, t_:  t] := \ 

      \ Module[{}, 
\ \ \ 
        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  =  {p'[t]  + w^2 x[t]  ==  0,  p[t] == x'[t]};

xfirstDiffSeries[eq1, {x, p}, {x[0] -> x0, p[0] -> p0}, 4, t];

Rearranged input equations with right side == 0 .

eq

Perturbative expansion of position:

zSer

Perturbative expansion of momentum:

eSer

List of unknown terms to be solved for:

vars

Solutions for these unknown terms:

sol

zSer and eSer with substitutions made:

zpSol

Example P1:  Perturbed harmonic oscillator.

eq1  =  
      x''[t] == -ω^2 x[t]  +  ϵ x[t]^2

firstOrderPert[eq1, x, {x0, 0}, ϵ, t]

How firstOrderPert works.

\(xfirstOrderPert[eq1_, x_, {x0_, v0_}, eps_:  eps, t_:  t] := 
\ \ 
      Module[{}, 
\ \ \ \ \ 
        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] == -ω^2 x[t]  +  ϵ x[t]^2

xfirstOrderPert[eq1, x, {x0, 0}, ϵ, t]

Expand the position variable:

x[t] /. xRule

Substitute this into modified input equation:

eq2

Group in powers of ϵ and set each == 0 .

eqEps

Unperturbed solution for position:

dsol1

Perturbed solution for position:

dsol2

Substitute these into perturbed expression for position:

x[t] /. xRule //. Join[dsol1, dsol2]