# Avoiding a singularity in the simulation of a spherica pendulum

I didn't know whether to put this here or in StackOverflow - so I open to answers just telling me to go there!

I am looking to simulate the motion of a spherical pendulum.

The Lagrangian is $$\frac{1}{2}(\dot{\theta}^2 + \sin^2 (\theta) \dot{\phi} ^2) + \cos (\theta)$$

where I have set $m=g=l=1$. $\theta$ and $\phi$ are the angular coordinates.

Since $\phi$ does not appear in the Lagrangian it is ignorable and we conclude that

$$\dot{\phi} = \frac{p_{\phi}}{\sin^2 \theta}$$

where $p_{\theta}$ is a constant associated with the $z-$component of the angular momentum.

Substituting this into an expression for the energy and rearranging we get, $$\dot{\theta} = \sqrt{2} \left(E+ \cos \theta + \frac{p_{\phi}}{2 \sin^2 \theta} \right)$$

where $E$ is the total energy.

I want to numerically integrate this system to find $\theta$ at any time $t$. Which integration method would be best for this? Would an RK45 avoid the division by zero when $\theta = n \pi$ ?

Thanks

The singular term is the centrifugal barrier and it is there to stop you from getting to $\theta=n\pi$. Once you fix an energy $E$ and an angular momentum $p_\phi$, the range of $\theta$ is restricted to the classically allowed region, between the two turning points.
Note, finally, that you will need to integrate with respect to $\theta$ and then invert that function - by interpolation if necessary - if you want to specify times $t$ and get the angle at those times. Depending on your application, though, you might be just fine with a list of times an angles $\{t_n,\theta_n\}$ which are evenly spaced in angle instead of time.