Using NDSolve to solve Equation of Motion in cylindrical coordinates

I have a set of coupled differential equations which represents the equation of motion of a particle in cylindrical coordinates with the following Hamiltonian: $$H=\frac{1}{2m} \left( p_r^2 + \frac{p_\phi^2}{r^2} + p_z^2 \right) + c \sqrt{\frac{r^2}{4 a^2} + \frac{z^2}{a^2} },$$ where a and c are some arbitrary constant.

From Hamilton's equations: $$\dot{q_i} = \frac{\partial H}{\partial p_i} \\ \dot{p_i} = - \frac{\partial H}{\partial q_i},$$

we get: $$\dot{r} = \frac{p_r}{m} \\ \dot{\phi} = \frac{p_\phi}{m r^2} \\ \dot{z} = \frac{p_z}{m} \\ \dot{p_r} = \frac{p_\phi^2}{m r^3} - \frac{c r}{4 a^2 \sqrt{\frac{r^2}{4 a^2} + \frac{z^2}{a^2} }} \\ \dot{p_\phi} = 0 \\ \dot{p_z} = - \frac{c z}{4 a^2 \sqrt{\frac{r^2}{4 a^2} + \frac{z^2}{a^2} }}$$

hamiltonianCylindrical[r_,phi_,z_,pr_,pPhi_,pz_,m_,a_,c_]:=1/(2 m) (pr^2+pPhi^2/r^2+pz^2) + c Sqrt[r^2/(4 a^2)+z^2/a^2]

cylindricalTrajectory[r0_, phi0_, z0_, pr0_, pPhi0_, pz0_, m_, a_, c_,steps_]:=NDSolve[
{
r'[t]==D[hamiltonianCylindrical[r[t],phi[t],z[t],pr[t],pPhi[t],pz[t],m,a,c],pr[t]],
phi'[t]==D[hamiltonianCylindrical[r[t],phi[t],z[t],pr[t],pPhi[t],pz[t],m,a,c],pPhi[t]],
z'[t]==D[hamiltonianCylindrical[r[t],phi[t],z[t],pr[t],pPhi[t],pz[t],m,a,c],pz[t]],
pr'[t]==-D[hamiltonianCylindrical[r[t],phi[t],z[t],pr[t],pPhi[t],pz[t],m,a,c],r[t]],
pPhi'[t]==-D[hamiltonianCylindrical[r[t],phi[t],z[t],pr[t],pPhi[t],pz[t],m,a,c],phi[t]],
pz'[t]==-D[hamiltonianCylindrical[r[t],phi[t],z[t],pr[t],pPhi[t],pz[t],m,a,c],z[t]],
r==r0,
phi==phi0,
z==z0,
pr==pr0,
pPhi==pPhi0,
pz==pz0
},
{r[t], phi[t], z[t], pr[t], pPhi[t], pz[t]},
{t, 0, steps},
WorkingPrecision -> 20,
MaxSteps -> \[Infinity]
]

For example, if I run the following:

Plot[
Evaluate[r[t]/.cylindricalTrajectory[1, 0, 2, 1, 0, 1, 1, 1, 1, 500]],
{t,0,500}
]

I can see that r becomes negative, which is not possible in cylindrical coordinates. Therefore, I wanted to know how I can solve the set of differential equation with r>=0.

Update

I noticed something weird is happening (which I have seen before):

Plot[
Evaluate[r[t] /. cylindricalTrajectory[1, 0, 2, 1, 1, 1, 1, 1, 1, 1000]],
{t, 0, 1000},
ImageSize -> Large
] At around t=500, r(t) behaves unexpectedly. However, if I only plot from t=0 to t=600, the "hiccups" at t=500 disappear.

Another question I have is how do I convert from Cartesian to Cylindrical, if I am correct, the following should be the transformation: $$r = \sqrt{x^2 + y^2} \\ \phi = arctan \left( \frac{y}{x} \right) \\ z=z \\ p_r = p_x cos \phi + p_y sin \phi \\ p_\phi = m (x p_y - y p_x) \\ p_z = p_z$$ The one I am not quite sure is $p_r$.

• I have no idea, but this may help visualising things: upt = 500; ct = cylindricalTrajectory[1, 0, 2, 1, 0, 1, 1, 1, 1, upt]; traj = Flatten[#, 1] &@ Table[{r[t], phi[t], z[t]} /. ct, {t, 0, upt, .1}]; trajcart = {#[]*Cos[#[]], #[]*Cos[#[]], #[]} & /@ traj; Manipulate[ Graphics3D[ { Arrow[{{0, 0, 0}, #}] &@trajcart[[i]], Point[trajcart[[1 ;; i]]] }, PlotRange -> {{-5, 5}, {-5, 5}, {-5, 5}} ], {i, 1, Length@trajcart, 1} ] – acl Jul 7 '14 at 20:16
• Maybe this helps with the glitch in the plot. – acl Jul 8 '14 at 12:04
• Ah, that fixed my plot. Thanks acl! – user29165 Jul 8 '14 at 12:16
• If you're not sure what the correct expression for the canonically conjugate momentum is, write down the Lagrangian in Cartesian coordinates, make the coordinate transformation, and then obtain the momentum from $p_r \equiv \partial{L}/\partial{\dot r}$. – Jens Jul 8 '14 at 16:12

The sign of the radial coordinate is of course supposed to be non-negative, and indeed your example shows that the sign flips to a negative value after some time. However, this is an artifact of the special case you considered in the example: the case of zero angular momentum, pPhi = 0.

If you give the angular momentum any arbitrarily small non-zero value, then the centrifugal barrier in the Hamiltonian will kick in and prevent r[t] from reaching the origin. Without that term, there is no way for Mathematica to tell from the equation that you mean to describe a polar coordinate system and not a Cartesian coordinate.

This is not a real problem because the negative-r solution for zero angular momentum is not strictly wrong. As the plot in acl's comment shows, the negative sign of r conspires with the unchanged polar angle phi to give the correct 3D trajectory. If one wanted to describe the same trajectory with a non-negative r then it would require flipping the polar angle from $\phi$ to $\phi+\pi$ at the instant the origin is crossed in real space, which however cannot be achieved by the Hamilton equations of motion because they require $\phi$ to remain constant if the conjugate momentum is zero. So what Mathematica outputs is really completely consistent.

Since I already have a solver for Hamilton's equations that I find easier to use than the code in the question, I'll repeat the same calculation with my code instead. The main function that solves the equations of motion is called hamiltonSolve, and it takes the expression for the Hamiltonian as its first argument. The second argument varList consists of pairs { {r, r0}, ... {z, z0}, ... {pz,pz0}, {t, tMax}}, where all the variable names are listed together with their initial values (quantities ending in 0). The last entry in the list is the name of the time variable, and its maximum value tMax, assuming the initial time is always t = 0.

hamiltonSolve[hamilton_, varList_] := Module[
{inputVariables, initialValues, phaseSpace, t, tMax, localH},
{t, tMax} = Last[varList];
{inputVariables, initialValues} = Transpose[Most[varList]];
phaseSpace = Through[inputVariables[t]];
localH = hamilton /. Thread[inputVariables -> phaseSpace];
inputVariables /. First[
NDSolve[
Map[D[#, t] &, phaseSpace] == Flatten[
Reverse[{-1, 1}*Partition[
Map[
D[localH, #] &,
phaseSpace
], Length[phaseSpace]/2
]
]
]],
Thread[(phaseSpace /. t -> 0) == initialValues
]]
,
inputVariables,
{t, 0, tMax},
MaxSteps -> 40000]
]
]

Now to apply this solver to your Hamiltonian, I repeat the definition from the question. Then I calculate the trajectory over a smaller time interval from 0 to 10, using the parameters in the question:

hamiltonianCylindrical[r_, phi_, z_, pr_, pPhi_, pz_, m_, a_, c_] :=
1/(2 m) (pr^2 + pPhi^2/r^2 + pz^2) + c Sqrt[r^2/(4 a^2) + z^2/a^2]

sol0 = With[{r0 = 1, phi0 = 0., z0 = 2, pr0 = 1, pPhi0 = 0., pz0 = 1,
m = 1, a = 1, c = 1},
hamiltonSolve[
hamiltonianCylindrical[r, phi, z, pr, pPhi, pz, m, a,
c], {{r, r0}, {phi, phi0}, {z, z0}, {pr, pr0}, {pPhi,
pPhi0}, {pz, pz0}, {t, 20}}]
];
Plot[sol0[][t], {t, 0, 10}, Axes -> True] Here I plotted the radial coordinate for zero angular momentum and indeed find the zero crossing. Now I simply change the angular momentum (which is a conserved quantity here) to pPhi = 0.001 and repeat the solution:

sol = With[{r0 = 1, phi0 = 0., z0 = 2, pr0 = 1, pPhi0 = 0.001,
pz0 = 1, m = 1, a = 1, c = 1},
hamiltonSolve[
hamiltonianCylindrical[r, phi, z, pr, pPhi, pz, m, a,
c], {{r, r0}, {phi, phi0}, {z, z0}, {pr, pr0}, {pPhi,
pPhi0}, {pz, pz0}, {t, 20}}]
];
Plot[sol[][t], {t, 0, 10}] The zero crossing has disappeared because the origin cannot be reached anymore. Other than that, the trajectory looks pretty much the same because the small angular momentum doesn't make any difference in the effective potential at larger radial distances.

Edit: polar angle

Following the suggestion by acl: To illustrate the singular character of the polar coordinate system when trying to describe the purely linear motion in the special case of zero angular momentum, you can look at the limit of small angular momentum and verify that the angle approaches a discontinuous jump. From the last solution, I pick the second function representing phi[t]:

Plot[sol[][t], {t, 0, 10}, PlotRangePadding -> .1] This shows that by enforcing $r \ge 0$ in this limit we will end up with a jump in $\phi$. But since Hamilton's equations tell us that $\phi = const$ at $p_{\phi} = 0$, Mathematica isn't doing anything wrong by obeying the latter equation.

• That's a great answer. Probably worth showing both $r$ and $\phi$, Plot[sol0[][t], {t, 0, 10}, Axes -> True] Plot[sol0[][t], {t, 0, 10}, Axes -> True], to explicitly show $\phi$ rapidly rotating when the radial component hits the centrifugal barrier (so the particle zooms by the origin). – acl Jul 7 '14 at 23:25
• @acl Good idea. I'll do it soon... – Jens Jul 8 '14 at 0:56
• Thanks Jens and acl, I agree Mathematica is solving the differential equations correctly. I have a couple more questions which I will update in my question above. – user29165 Jul 8 '14 at 11:12