The goal of the present paper can then be formulated as follows: find for each $\epsilon >0$ and each real $s,p$ $(p\geq 1)$ an approximative solution $q^2_{s,p}$ to such that $$\|d^2-q_{s,p}^2\|_{s,p}\leq \epsilon .$$ We shall call $q_{s,p}^2$ an $H^{s,p}$ approximation to $d^ 2$ for reasons which will become apparent later. Let us motivate this ambitious task by looking at a specific application. There are a lot of applications for computations of the length of a geodesic, where applications to computations in general relativity are only one domain. Another important example is the leading term of the expansion of the fundamental solution of linear parabolic solutions (with variable coefficients). Varadhan showed that the fundamental solution of the diffusion equation $$\label{PPDE} 
\begin{array}{l}
\frac{\partial u}{\partial t}=\frac{1}{2}\sum_{i,j}a_{ij}\frac{\partial^2 u}{\partial x_i\partial x_j}+
\sum_i b_i\frac{\partial u}{\partial x_i},
\end{array}$$ (where the diffusion coefficients $a_{ij}$ and the first order coefficients $b_i$ in depend on the spatial variable $x$ only) is connected to the length $d$ of the geodesic with respect to the line element $ds^2=\sum_{ij}a^{ij}dx_idx_j$ ($a^{ij}$ being the inverse of $a_{ij}$) via the relation $$d^2(x,y)=\lim_{t\downarrow 0}t\ln p(t,x,y).$$

Solving equation we can assume that the matrix-valued function $x\rightarrow (a_{ij}(x))$ is symmetric, i.e. $a_{ij}(x)=a_{ji}(x)$ for all $1\leq i,j\leq n$. This is 