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.2  Orbit Problems

Problem 4:  Numerical Solution for Orbits with Central Forces

Clear["Global`*`"];

Part (a)

eq1=u''[ϕ]==-u[ϕ] - 
            ((m f[(1/u)[ϕ]])/((el^2) u[ϕ]^2));

forceRule={\(V -> \((\(-\(k\/#1\)\)&)\)\),f[(1/u)[ϕ]]->-V'[r],r -> (1/u)[ϕ]};

eq2 = eq1 //. forceRule

values = {k -> 1, el -> 1, m -> 1};

initial1={u[0] == 1.6,u'[0]==0};

ndsol1 = 
      Flatten[NDSolve[Flatten[{eq2 /. values, initial1}], 
          u[ϕ], {ϕ, 0, 2 π}]];

u1[ϕ_] = u[ϕ] //. ndsol1; u1[1.1]

PolarPlot[Evaluate[(1/u1)[ϕ]], {ϕ, 0, 2 π}];

initial2={u[0] == 3,u'[0]==0};

ndsol2 = 
      Flatten[NDSolve[Flatten[{eq2 /. values, initial2}], 
          u[ϕ], {ϕ, 0, 2 π}]];

u2[ϕ_] = u[ϕ] //. ndsol2; u2[1.1]

PolarPlot[Evaluate[(1/u2)[ϕ]], {ϕ, 0, 2 π}];

Clear[numOrbit]; numOrbit[V_, r_, r0_:  1, dr0_:  0, el_:  1, m_:  1, 
        angle_:  6 π]:=Module[{u, ϕ, uforce, eq1, initial, eqs, ndsol},uforce = -(∂_r) V /. r -> (1/u)[ϕ];eq1=ExpandAll[u''[ϕ]==-u[ϕ] - 
                    ((m uforce)/((el^2) u[ϕ]^2))];initial={u[0] == (1/r0),u'[0]==-(dr0/(r0^2))};eqs = Join[{eq1}, initial];ndsol = Flatten[NDSolve[eqs, u[ϕ], {ϕ, 0, angle}]];PolarPlot[
              Evaluate[(1/u)[ϕ] /. ndsol], {ϕ, 0, angle}]]

numOrbit[-(1/r), r, 0.6];

numOrbit[-(1/(r^1.1)), r, 0.6];

numOrbit[-(1/(r^0.9)), r, 0.6];

numOrbit[-(1/r) - (0.02/(r^3)), r, 0.6];

numOrbit[-(1/r) + (0.02/(r^3)), r, 0.6];