want to compute. If we knew, then computing the length would be a rather trivial task. Even the derivatives of the length functional would be better computed from the explicit geodesic. However, as we mentioned the nonlinear ordinary differential equation describing the geodesic is harder to solve in general than the eikonal equation. Second, we shall see from an different point of view why the boundary condition leads to uniqueness of solutions $(x,y)\rightarrow d^2(x,y)$ of the system , , and . We have

The local representation $d^2$ satisfying the equations , , together with the boundary condition is of the form $$\label{loc}
\begin{array}{ll}
d^2(x,y)&=\sum_{ij}a^{ij}(y)\Delta x^i\Delta x^j + \sum_{|\alpha < M }\frac{d^2_{\alpha}(y)}{\alpha!}\Delta x^{\alpha}\\
\\
&+\sum_{|\gamma|=M}\int_0^1(1-\theta)^{M-1}\frac{\Delta x^{\gamma}}{\gamma !}\partial^{\gamma}d^2(y+\theta \Delta x,y)d\theta.
\end{array}$$ The coefficients $d_{\alpha}(y)$ are uniquely determined by a recursion obtained from the equations , . In coordinates with second order normal form, i.e. where $d^2$ is $\sum_{ij}\lambda_i(y)\Delta x^i\Delta x^j$ with $\lambda_i(y) , 1\leq i \leq n$ is the spectrum of $(a^{ij}(y))$, the multiindex recursion is $$\label{rec}
\begin{array}{ll}
 d^2_{\beta}(y)=&\frac{1}{\left(1-\sum_i\beta_i\right)}{\Bigg (} \sum_i  \left(  \lambda^i_0\right)^2\frac{\lambda_i^{\beta\dot -2_i}}{(\beta\dot -2_i)!}1_{\left\lbrace \beta_i\geq 2\right\rbrace }\\
\\
&+\sum_i \sum_{|\alpha|\geq 1,|\gamma | \geq 3, \alpha+\gamma =\beta}\frac{\lambda_i^{\alpha}}{\alpha!}  \lambda^i_0 d^2_{\gamma}(y)\gamma_i\\
\\
&+ \sum_i  \sum_{\alpha\geq 0,|\delta|\geq 3,|\gamma| \geq 3,\alpha+\gamma+\delta\dot-2_i=\beta}\frac{\lambda_i^{\alpha}}{\alpha!}  \delta_i \gamma_i d^2_{\delta}(y)d^2_{\gamma}(y){\Bigg )}.
\end{array}$$ 