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):


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$):

enter image description here

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:

enter image description here

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

  • $\begingroup$ 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. $\endgroup$
    – user10851
    Aug 27, 2013 at 18:55
  • 1
    $\begingroup$ @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. $\endgroup$
    – Ruslan
    Aug 27, 2013 at 19:21

1 Answer 1


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}$$

With this corrected, the result is as expected:

enter image description here


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.