# Numerical integration of differential equation on the surface of a sphere

I'm trying to simulate the motion of a particle (a position vector) that is constrained to living on the surface of a unit sphere. Each time-step, $$\Delta t$$, the particle moves in some direction on the sphere's surface according to some representation of velocity/ acceleration/ etc. Think Euler method in spherical coords. I feel like this should be a straight-forward task, but I'm suffering a major mind block. In fact, I'm not even sure if I'm posing the question correctly.

To state it differently, given two previous positions on a unit sphere, $$\vec{r}_{t-1}$$ and $$\vec{r}_{t}$$, integrate forward to find the next position $$\vec{r}_{t+1}$$, assuming linear dynamics (i.e., the particle is moving in a straight "line" on the surface of the sphere).

To test if the model is implemented correctly, I should be able to select an arbitrary point on the sphere, give it a constant velocity in some arbitrary direction, and it should circle around the sphere and return to its starting position.

What's the optimal way to do this? Should the particle's velocity be in terms of $$\theta$$ (colatitude) and $$\phi$$ (azimuth), or something else? Thanks for any pointers!

• If your system can be described with a Hamiltonian, a symplectic integrator will work better on long times, and has an easy geometric description that will allow you to work on any coordinate system. – mlainz Jun 3 at 18:43
• @R.Taylor Search for "spherical pendulum". You can derive a system of differential equations in cartesian coordinates that will describe the motion of a particle in a spherical surface. The "test" you propose is only valid if you discard external forces like gravity, but the solutions are far more interesting if you do not. – PierreCarre Jun 4 at 21:39

If you denote by $$X(t)$$ the position of the particle at time $$t$$, the set of differential equations that describe the motion is
$$\ddot X=\frac 1m \left( F - \dfrac{m \dot X^T H \dot X + \ \nabla \Phi^T F}{\|\nabla \Phi \|^2} \nabla \Phi\right), \quad t >0,$$
where $$\Phi(x) = x_1^2+x_2^2+x_3^2-1$$, $$\nabla \Phi$$ is its spatial gradient and $$H$$ its Hessian matrix.