Here I describe a procedure for computing random walks on surfaces (that is, manifolds of dimension two). The procedure can be easily generalised to higher dimensions, but surfaces are more commonly encountered in practice.
Let $\vec x(u,v) = (x(u,v),y(u,v),z(u,v))$ represent a parameterised surface, and $(u_0,v_0)$ the coordinates of the initial position of a diffusing particle. We illustrate with the example of a torus of major radius $R$ and minor raidus $r$, which can be parameterised by \[ \vec x = [\cos{u}(R+r\cos{v}),\sin{u}(R+r\cos{v}),r\sin{v}]. \] Of course, we cannot simply perform a random walk in the coordinate space $(u,v)$. Not only would this depend on the choice of parameterisation, but in general there exists no parameterisation for which a random walk in the coordinate space coincides with a random walk on the actual surface.
Instead, in order to perform a random walk on a surface, we need to find a pair of vectors $e_1(u,v), e_2(u,v)$ in the coordinate space which represent orthonormal directions on the surface. To do this, we begin by computing the Riemannian metric, which is defined by the following matrix: \[{g_{ij}} = \left[ {\begin{array}{*{20}{c}} {{g_{11}}}&{{g_{12}}}\\ {{g_{21}}}&{{g_{22}}} \end{array}} \right] = \left[ {\begin{array}{*{20}{c}} {\frac{{\partial \vec x}}{{\partial u}} \cdot \frac{{\partial \vec x}}{{\partial u}}}&{\frac{{\partial \vec x}}{{\partial u}} \cdot \frac{{\partial \vec x}}{{\partial v}}}\\ {\frac{{\partial \vec x}}{{\partial v}} \cdot \frac{{\partial \vec x}}{{\partial u}}}&{\frac{{\partial \vec x}}{{\partial v}} \cdot \frac{{\partial \vec x}}{{\partial v}}} \end{array}} \right] = \left[ {\begin{array}{*{20}{c}} {{{(R + r\cos{v})}^2}}&0\\ 0&{{r^2}} \end{array}} \right].\] Using the Gramm-Schmidt process, we can take our vectors to be \begin{align*} e_1 & = \frac{1}{\sqrt{g_{11}}}(1,0),\\ e_2 & = \frac{1}{\sqrt{g}\sqrt{g_{11}}}(-g_{21},g_{11}), \end{align*} where $g$ is the determinant of the above matrix. In the specific case of the torus, this works out to be \begin{align*} e_1 & = ((R+r\cos{v})^{-1},0), \\ e_2 & = (0,r^{-1}). \end{align*} Now, all we have to do to choose at random a vector $e_0$ from the four vectors ${\pm e_1, \pm e_2}$, and the new coordinates of the particle will be given by \[ (u_1,v_1) = (u_0,v_0) + \delta e_0, \] where $\delta$ is a small number representing the step size.
