How would one discretize the dynamical equations for Kerr spacetime?

Question

Using the Hamiltonian for a test particle in Kerr spacetime, we arrive at the following equations for generalized position and momenta (in natural units, $G = c = M = 1$):

\begin{align} \dot{r} &= \frac{\Delta}{\Sigma}p_r\\ \dot{\theta} &= \frac{1}{\Sigma}p_\theta\\ \dot{\phi} &= \frac{1}{\Sigma\Delta}\left(2arE + (\Sigma - 2r)L\csc^2\theta\right)\\ \dot{p}_r &= \frac{1}{\Sigma\Delta}\left[(r - 1)((r^2 + a^2)\mu - \kappa + r\Delta\mu + 2r(r^2 + a^2)E^2 - 2aEL\right]-\frac{2p_r^2(r-1)}{\Sigma}\\ \dot{p}_\theta &= \frac{\sin\theta\cos\theta}{\Sigma}\left[L^2\csc^4\theta - a^2(E^2 + \mu)\right]\\ \dot{p}_\phi &= 0 \end{align} where, $\kappa = p_\theta^2 + L^2\csc^2\theta + a^2 - (E^2\sin^2\theta + \mu)$, $\Sigma = r^2 + a^2\cos^2\theta$, $\Delta = r^2 - 2r + a^2$; and $\mu$, $E = -p_t$ and $L = p_\phi$ are the test particle mass, energy and angular momentum, respectively. Also, $\mu = -1$ for massive particles, while for massless particles, $\mu = 0$, and the derivatives are with respect to an affine parameter, say $\lambda$.

I want to implement a symplectic leapfrog/verlet solver to solve this system. But I don't understand, how this set of equations should be discretized. Since $p_i$ are the momenta and not velocities, does the usual Taylor series approach make sense here? If not, could someone help me understand, how a verlet/leapfrog scheme can be used here?