$N$ goes to infinity in such a way that the mesh size of the set of interpolation points $h$ goes to zero. Note that $q_{0,p}$ denotes the squareroot of $q_{0,p}^2$. We call an approximation $q^2_{s,p}$ an $H^{s,p}$-approximation if it approximates not only $d^2$ in the $L^p$ sense but can be plugged in into all the derivatives of of order $m$ (i.e. multivariate derivatives $\alpha$ for $|\alpha|\leq  m$ of the eikonal equation) such that in $$\frac{\partial^{\alpha}}{\partial x^{\alpha}}\left(\frac{1}{4}\sum_{ij}a_{ij}(x)\frac{\partial q^2_{0,p}}{\partial x_i}\frac{\partial q^2_{0,p}}{\partial x_j} \right)-\frac{\partial^{\alpha}}{\partial x^{\alpha}}d^2(x,y)=:r_{\alpha ,p}$$ the right side staisfies $r_{0,p}\in O(h^{3+m})$ locally implies that $$\|d(x,y)-q_{0,p}\|_{H^{s,p}(\Omega)}$$ converges to zero as the number of interpolation points $N$ goes to infinity in such a way that the mesh size of the set of interpolation points $h$ goes to zero. Accordingly, we call such $q_{0,p}^2$ ($q_{s,p}^2$) an $L^p$- ($H^{s,p}$) approximation of the boundary value problem . In the next subsection we construct a $L^p$-approximation and refine the construction in the following subsection in order to construct $H^{s,p}$-approximations.

Polynomial interpolation of eikonal equation in $L^ p$ sense
------------------------------------------------------------

We may write the eikonal equation $$\label{eik0}
d^2(x,y)=\frac{1}{4}\sum_{ij}a_{ij}d^ 2_{x_i}d^ 2_{x_j}=\frac{1}{4}\left\langle \nabla d^2, A\nabla d^2\right\rangle.$$ Assume that $A=(a_{ij})$ is constant. The solution with the boundary condition $d^2(x,y)=0$ iff $x=y$ is 