# Two interacting particles on sphere drift to sphere poles

Suppose we have two particles which can move on sphere of radius $r$, and they attract to each other so that their potential energy is $g(d)=ad$ where $d$ is distance between them. I've found Lagrangian, it looks like this (in spherical coordinates):

$$L=\frac{r^2}2\left(m_1\left(\dot\theta_1^2+\dot\varphi_1^2\sin^2\theta_1\right)+m_2\left(\dot\theta_2^2+\dot\varphi_2^2\sin^2\theta_2\right)\right)-g\left(l_{arc}\left(\theta_1,\varphi_1,\theta_2,\varphi_2\right)\right),$$

where $l_{arc}$ is arc length between two points on sphere: $$l_{arc}=2r\arcsin\left(\frac1{\sqrt2}\sqrt{1-\cos\theta_1\cos\theta_2-\cos\left(\varphi_1-\varphi_2\right)\sin\theta_1\sin\theta_2}\right)$$

So, equations of motion for particle $i$ look like:

$$\frac{m_i r^2}2 \frac{\text{d}}{\text{d}t}\begin{pmatrix}2\dot\theta_i\\ 2\dot\varphi_i\sin^2\theta_i\end{pmatrix}=-\begin{pmatrix}\frac{\partial g(l_{arc})}{\partial \theta_i}\\ \frac{\partial g(l_{arc})}{\partial \varphi_i}\end{pmatrix}+\begin{pmatrix}\sin\left(2\theta_i\right)\dot\varphi_i^2\\ 0\end{pmatrix}$$

I solve a Cauchy problem, so here're initial conditions: $$\theta_1(0)=\frac\pi2+10^{-4}\\ \theta_2(0)=1.05\cdot\frac\pi2\\ \varphi_1(0)=-1.5\\ \phi_2(0)=-1.45\\ \dot\theta_1(0)=0.003\\ \dot\theta_2(0)=0.003\\ \dot\varphi_1(0)=-0.01\\ \dot\phi_2(0)=-0.01$$

Now as I solve this problem, it appears that the system drifts to the pole of the sphere whenever the particles interact (here $r=100,\; m_1=m_2=1,\; a=1$):

I've tried multiple methods of solving this problem in Mathematica using NDSolve with big range of steps, but still the trend is the same, it seems to not depend on method of solution. So, I seem to have made a mistake somewhere in formulating the problem.

If instead of attracting the particles repulse, they again appear attracted by poles:

Is my derivation correct - i.e. finding Lagrangian and deriving equations of motion? What could be the reason for such strange error?

-
Singling out the coordinate poles when the physical situation says they're not special is a sign of a numerical problem, not a physical one. Mathematica is great for symbolic manipulation, but it is really not good at all at numerical stability issues. There's a whole field dedicated to these sorts of things precisely because you often can't throw a prepackaged, blackbox routine at a simulation and expect the results to be physically meaningful. –  Chris White Aug 27 '13 at 18:55
@ChrisWhite yeah, your statement makes sense, but the problem may be wrongly formulated, and this could also lead to poles becoming special. In fact I've now found the mistake, it was just omitted factor of $\frac{mr^2}2$ in the rightmost vector in equations of motion. I'm gonna answer my own question. –  Ruslan Aug 27 '13 at 19:21

My mistake was omitting factor of $\frac{mr^2}2$ in equations of motion near rightmost vector. They should instead look like this:
$$\frac{m_i r^2}2 \frac{\text{d}}{\text{d}t}\begin{pmatrix}2\dot\theta_i\\ 2\dot\varphi_i\sin^2\theta_i\end{pmatrix}=-\begin{pmatrix}\frac{\partial g(l_{arc})}{\partial \theta_i}\\ \frac{\partial g(l_{arc})}{\partial \varphi_i}\end{pmatrix}+\frac{m_i r^2}2\begin{pmatrix}\sin\left(2\theta_i\right)\dot\varphi_i^2\\ 0\end{pmatrix}$$