space, and it is well known that for any compact Riemannian manifold any two points $x,y\in M$ can be connected by a geodesic whose length is $d_M(x,y)$. If $\nabla$ denotes the Levi-Civita connection, then a geodesic $\gamma$ is characterized by the equation $$\nabla_{\dot\gamma}\dot\gamma=0,$$ which becomes (in terms of the coordinates of the values of the curve $\gamma$) $$\label{geo1}
    \frac{d^2x^\lambda }{dt^2} + \Gamma^{\lambda}_{~\mu \nu }\frac{dx^\mu }{dt}\frac{dx^\nu }{dt} = 0\ ,$$ where the well-known Christoffel symbols are $$\label{geo2}
    \Gamma^\kappa_{\; \mu \nu}=\frac{1}{2}g^{\kappa \rho} \left( \partial_\mu g_{\nu \rho}+\partial_\nu g_{\mu \rho}-\partial_\rho g_{\mu \nu} \right).$$ This is an $n$-dimensional nonlinear ordinary differential equation with values in ${\mathbb R}^n$ which is difficult to compute numerically in general (note the quadratic terms). For computing the length of a geodesic it is easier to compute the solution of a eikonal equation of the form $$\label{eik}
d^2=\frac{1}{4}\sum_{ij}a_{ij}(x)d^2_{x_i}d^2_{x_j}$$ (boundary conditions considered later), where $x\rightarrow a_{ij}(x)$ are functions such that at each $x\in{\mathbb R}^n$ the matrix $(a_{ij}(x))$ is the inverse of the positive matrix $(g_{ij}(x))$ at each point $x$. Here $f_{x_i}:=\frac{\partial f}{\partial x_i}$ denotes the derivative of $f$ with respect to the variable $x_i$. In general we shall write $\partial^{\alpha}f$, $\partial^{\alpha}_xf$ or $\frac{\partial}{\partial x^{\alpha}}f$ for the multivariate derivative with multiindex $\alpha=(\alpha_1,\cdots ,\alpha_n)$. 