This notebook contains additional material to the Mathematica GuideBooks  (Springer-Verlag, New York, 2004) http://www.MathematicaGuideBooks.org . It uses the stylesheet of the Mathematica GuideBooks, which is embedded.

Symbolics Volume, Chapter 1

Given a 1D holonomic mechanical system with a conserved quantity K(x(t), v(t)), (v(t)=x'(t)), the Kobussen–Leubner–Lopez [✶] formula allows for the construction of a Lagrangian L(x,v) through L(x,v)=v ∫K(x,v)/(v^2)dv. (Given a Lagrangian, unique equations of motion follow, but given the equations of motion, the Lagrangian is nonunique [✶].) For a one-dimensional harmonic oscillator, the energy ε=m (v^2)/2+m (ω^2)(x^2)/2 is a conserved quantity. Using the Kobussen–Leubner–Lopez formula for powers of the energy (ε^n) (n∈(SuperPlus[ℕ])) yields the Lagrangian [✶]

(Subscript[L, n])(x,v)=(m/2)(∑_(k=0)^n)(1/(2k-1))(({{n}, {k}}))((((ω^2)(x^2)))^(n-k))(v^(2k)).

The solution of the Euler–Lagrange equations of (Subscript[L, n])(x,v) are the well known harmonic oscillations of the one-dimensional harmonic oscillator. But the phase space trajectories derived from the Hamiltonian (Subscript[H, n])(x,p), where (Subscript[H, n])(x,p)=((ε(x, v(x,p)))^n), and p=∂L/∂v are for n>1 different curves than the well-known circles of the one-dimensional harmonic oscillator [✶]. For n>1, calculate and visualize some of the phase space trajectories following from (Subscript[H, n])(x,p) via x'(t)=∂H(x(t),p(t))/∂p(t) and p'(t)=-∂H(x(t),p(t))/∂x(t).

The function implicitHamiltonianEquation calculates the implicit equation for the Hamiltonian. (We set m=ω=1 in the following.)

implicitHamiltonianEquation[n_, H_, {x_, p_}] := 
With[{L = 1/2 Sum[Binomial[n, k] (x^2)^(n - k) v^(2k)/(2k - 1), {k, 0, n}]},
     (* eliminate v *)
     Resultant[H - (1/2)^n (v^2 + x^2)^n, p - D[L, v], v] // 
                                             Together // Factor // Numerator]

Here are the first three polynomials (Subscript[𝓅, n])(x,p,H).

Table[implicitHamiltonianEquation[n, H, {x, p}], {n, 3}]

As described, the equations of motions are obtained from implicit differentiation.

HamiltonianEquationsOfMotion[n_, H_, {x_, p_}, t_] :=
With[{𝓅 = implicitHamiltonianEquation[n, H[x[t], p[t]], {x[t], p[t]}]},
     (* differentiate implicitly, then solve for derivatives of H *)
     {x'[t] == +Solve[D[𝓅, p[t]] == 0, D[H[x[t], p[t]], p[t]]][[1, 1, 2]],
      p'[t] == -Solve[D[𝓅, x[t]] == 0, D[H[x[t], p[t]], x[t]]][[1, 1, 2]]}]

For n=1, we recover the well-known equations of motion for the harmonic oscillator.

HamiltonianEquationsOfMotion[1, H, {x, p}, t]

For n=2, we obtain considerably more complicated equations. They contain the Hamiltonian H, a conserved quantity.

HamiltonianEquationsOfMotion[2, H, {x, p}, t]

Now, we must deal with the initial conditions for x and p. We will prescribe values for x(0) and p(0). This fixes the value of H, but not necessary uniquely, because H appears in (Subscript[p, n]) nonlinearly. We use high-precision arithmetic to calculate the starting values and to solve the differential equaion to obtain numerically reliable solutions.

HValues[n_, {x0_, p0_}, wp_: 30] :=
Module[{H, NHs, x, p, t, x0E = Rationalize[x0, 0], 
        p0E = Rationalize[p0, 0],dens}, 
       (* numerically solve for H *)
       NHs = N[Solve[implicitHamiltonianEquation[n, H, {x0E, p0E}] == 0, H], wp];
       (* eliminate potential zeros in the denominators of the derivatives of H *)
       dens = Denominator[Last /@ HamiltonianEquationsOfMotion[n, H, {x, p}, t]] //.
                   ξ_[t] :> ξ /. _H :> H;      
       (* esort in such a manner that real solutiosn come first *)
       Sort[#[[1, 2]]& /@ Select[NHs, FreeQ[dens /. # /. {x -> x0E, p -> p0E},
                                       _?(# == 0&)]&], Abs[Im[#1]] < Abs[Im[#2]]&]]

Generically, for real x(0), p(0), we have one real value for H. Here is an example.

HValues[4, {2, 2}] // N

Now, we are ready to solve the equations of motion numerically.

Clear[nDSolveHO];
nDSolveHO[n_, {x_, p_}, {x0_, p0_}, T_, k_:All, wp_:20] :=
Module[{hs},
       (* starting values *)
       hs = HValues[n, {x0, p0}, wp + 10];
       (* solve equations of motion for all starting values *)
       NDSolve[Join[HamiltonianEquationsOfMotion[n, H, {x, p}, t] /. H[__] -> #,
                    {x[0] == x0, p[0] == p0}],
               {x, p}, {t, 0, T}, WorkingPrecision -> wp, MaxSteps -> 10^5,
               PrecisionGoal -> 10, AccuracyGoal -> 10]& /@ Take[hs, k]]

To visualize the phase space trajectories, we use the function phaseSpaceTrajectories. It generates graphics of the real and imaginary parts of the phase space trajectories.

phaseSpaceTrajectories[n_, {x0_, p0_}, T_, k_:All, wp_:20] := 
Module[{ndsols = nDSolveHO[n, {x, p}, {x0, p0}, T, k, wp]},
        λ = Length[ndsols]; 
       GraphicsArray[Function[reIm, 
           (* plot real and imaginary parts of the parametrized trajectories *)
           ParametricPlot[Evaluate[reIm @ {x[t], p[t]} /. ndsols],
                          {t, 0, T}, PlotRange -> All, 
                          PlotPoints -> Round[T 100], AspectRatio -> 1,
                          DisplayFunction -> Identity, PlotLabel -> reIm,
                          Frame -> True, Axes -> False,
                          PlotStyle -> Table[Hue[0.8 j/λ], 
                                                 {j, 0, λ - 1}]]] /@ {Re, Im}]]

For n=1, we have just one possible phase space trajectory, a circle (an ellipse).

Show @ phaseSpaceTrajectories[1, {2, 2}, 2Pi]

For n=2, we have three trajectories, one real one, and two coinciding complex one. The real one is displayed in red. The period is now doubled. The complex trajectory is more complicated.

Show @ phaseSpaceTrajectories[2, {2, 2}, 4Pi]

For n=2, we have five trajectories, one real one, and two pairs of coinciding complex one. The real one is again displayed in red. The period is four times the one from the n=1 case.  Again, the complex trajectories are more complicated.

Show @ phaseSpaceTrajectories[3, {2, 2}, 8Pi]

For n=3, we have seven trajectories, one real one, and three pairs of coinciding complex one. The real one is again displayed in red. The period is eight times the one from the n=1 case. And again, the complex trajectories are more complicated.

Show @ phaseSpaceTrajectories[4, {2, 2}, 16Pi]

We end here with n=4. Higher order complex trajectories become more and more complicated.

We end with displaying the phase space trajectories corresponding to the real trajectories for n=1,2,3 in one graphic. The different shape and size of the trajectories is clearly visible.

Show[Table[phaseSpaceTrajectories[n, {1/100, 1}, 2^n 2Pi, 1][[1, 1]] /. 
                          (* recolor to red (n = 1) to blue (n = 4) *) 
                                  Hue[0] :> Hue[(n - 1)/4], {n, 1, 4}],
     DisplayFunction -> $DisplayFunction, PlotRange -> All]

We could now continue and use complex starting values for x(0) and p(0) and solve the equations of motion in the complex t-plane to obtain still more complicated trajectories.