numerically integrating a trajectory in polar coordinates So I've reduced my problem to not being sure how to integrate a trajectory in polar coordinates.  Suppose I have a free particle and I express its Hamiltonian thus:
$H =\eta_{ij}P^iP^j,$
where $\eta_{ij}$ is the flat-space metric in 2D polar coordinates:
$
\eta_{ij} = \begin{pmatrix}
    1       & 0\\0&r^2 
\end{pmatrix}$.
And the $P^i$ vector contains the momenta: $P_i = (P_r, P_\theta)$
This gives a Hamiltonian of 
$H = P_r^2 + r^2P_\theta^2$.
(I know this isn't the most straightforward way to do this, but it generalizes to my problem in a useful way– I need to determine a first-order system for my equations of motion via the Hamiltonian given a metric).  Getting the equations of motion gives
$\dot{r} = 2P_r, \dot{P_r} = -2P_\theta^2 r,$
$\dot{\theta} = 2P_\theta r^2, \dot{P_\theta} = 0$.
Now, when I plug this into my integrator (a simple RK4 implementation where the initial condition is $(r,P_r,\theta,P_\theta)$), I get bizarre plots (this is the particle's trajectory, the axes are $x$ and $y$ (i.e. $r\sin(\theta), r\cos(\theta)$):

This "orbit" is not numerical error.  Increasing the time step by two orders of magnitude just makes it cleaner:

(interestingly, the initial condition for these plots was $r = 10, P_r = 1, \theta = \pi/3, P_\theta = 1$)
What's going on here? It should be a straight line, right?
 A: Your Hamiltonian doesn't mean what you think it means, because the canonical momenta are the covariant components of the momentum, not the contravariant ones.
If we start with the Lagrangian
$$L = g_{ij} \dot{q}^i \dot{q}^j$$
we can straightforwardly find that the canonical momenta are
$$p_i = 2 g_{ij} \dot{q}^j.$$
You can see by the index position on the RHS that it makes sense for $p_i$ to have a downstairs index, and that (modulo a normalization) it's nothing more than the covariant version of the velocity. After some index gymnastics we arrive at the Hamiltonian
$$H = g^{ij}p_i p_j.$$
Explicitly for your example, we have $L = \dot{r}^2 + r^2 \dot{\theta}^2$ and $H = p_r^2 + p_\theta^2/r^2$. If you use this Hamiltonian you should get the right trajectories.
There is, however, a subtlety here: My Hamiltonian and yours are equal! After all, index position doesn't matter in a contraction: $g^{ij}p_i p_j$ and $g_{ij} p^i p^j$ are the same. The problem, though, is that the $p^i$ are not canonical variables, so Hamilton's equations don't apply. To use the Hamiltonian formalism you need to use the covariant momentum.
