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?