Take the 2-minute tour ×
Physics Stack Exchange is a question and answer site for active researchers, academics and students of physics. It's 100% free, no registration required.

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

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?

share|improve this question
    
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
1  
@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

1 Answer 1

up vote 1 down vote accepted

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

share|improve this answer

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

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