1 Introduction
Many applied problems lead to the need of solving a BVP for the eikonal equation. First of all, this nonlinear partial differential equation is used to simulate wave propagation in the approximation of geometric optics
Born ; Kravtsov . In computational fluid dynamics, image processing and computer graphics (see, for example, sethian1999level ; Gilles ), the solution of BVPs for the eikonal equation is associated with calculating the nearest distance to boundaries of a computational domain.The eikonal equation is a typical example of steadystate Hamilton–Jacobi equations. The issues of the existence and uniqueness of the solution for boundary value problems for such equations are considered, e.g., in kruvzkov1975generalized ; lions1982generalized . To solve numerically BVPs for the eikonal equation, the standard approaches are used, which are based on using difference methods on rectangular grids or finiteelement/finitevolume approximations on general irregular grids. In this approach, the main attention is paid to problems of nonlinearity.
A boundary value problem is formulated in the following way. The function is defined as the solution of the equation
(1) 
in a domain with the specified boundary conditions
(2) 
Computational algorithms for solving BVPs for the eikonal equation can be divided into two classes.
Marching methods (the first class of algorithms) are the most widely used. They are based on the hyperbolic nature of the eikonal equation. In this case, the desired solution of the problem (1), (2) is obtained by successive moving into the interior of the domain from its boundary, using, for instance, firstorder upwind finite differences Tsitsiklis ; sethian1999fast . Among other popular methods, we should mention, first of all, the fast sweeping method tsai2003fast ; Zhao , which uses a Gauss–Seidelstyle update strategy to progress across the domain. Recently (see, for example, jeong2008fast ), a fast iterative method for eikonal equations is actively developed using triangular fu2011fast and tetrahedral fu2013fast grids. Other modern variants of the fast marching method, which are adapted, in particular, to modern computing systems of parallel architecture, have been studied and compared, e.g., in gomez2015fast .
The second class of algorithms is associated with a transition from (1), (2) to a linear or nonlinear BVP for an elliptic equation belyaev2015variational . Instead of equation (1), we (see li2005level ) minimize the functional
It is possible to solve the BVP for the Euler–Lagrange equation for this functional, which has the form
In fares2002differential , the computational algorithm is based on solving the nonlinear boundary value problem for . The solution of the problem (1), (2) can be related to the solution of the homogeneous Dirichlet problem for Laplacian:
In this case (see, e.g., bhattacharya1989limits ; kawohl1990family ), we have
Thus, to find the solution of the problem (1), (2), we need to solve the nonlinear BVPs.
In our study, we focus on solving auxiliary boundary value problems for linear equations. This approach (see gurumoorthy2009schrodinger ; sethi2012schrodinger )is based on a connection between the nonlinear Hamilton–Jacobi equation and the linear Schrodinger equation. Let be the solution of the boundary value problem
(3) 
(4) 
Then, for , we have as . A similar approach, where the auxiliary functions are associated with the solution of the unsteady heat equation, is considered in the paper crane2013geodesics .
In the present paper, we consider the eikonal equation in an anisotropic medium that is a more general variant in comparison with (1). Using the transformation , the corresponding BVP of type (3), (4) is formulated for the new unknown quantity. In our case, and so, we have a singularly perturbed BVP for the diffusion–reaction equation roos2008robust ; miller2012fitted . The numerical solution is based on using standard Lagrangian finite elements brenner2008mathematical ; larson2013finite . The main attention is paid to the monotonicity of the approximate solution for the auxiliary problem.
The paper is organized as follows. A boundary value problem for the eikonal equation in an anisotropic medium is formulated in Section 2. Its approximate solution is based on a transition to a singularly perturbed diffusion–convection equation. In Section 3, an approximation in space is constructed using Lagrangian finite elements and the main features of the problem solution are discussed. Numerical experiments on the accuracy of the approximate solution are presented in Section 5 for model twodimensional problems. The results of the work are summarized in Section 5.
2 Transformation of BVP for the eikonal equation in an anisotropic medium
In a bounded polygonal domain , with the Lipschitz continuous boundary , we search the solution of the BVP for the eikonal equation
(5) 
Define the operator as
(6) 
with the coefficients . The equation (5) is supplemented with the homogeneous Dirichlet boundary condition
(7) 
The basic problems of numerical solving the boundary value problem (1)–(3) result from the nonlinearity of the equation (see the operator ).
Similarly to gurumoorthy2009schrodinger ; sethi2012schrodinger , we introduce the transformation
(8) 
with a numerical parameter . This type of transformation is widely used in studying differential equations with quadratic nonlinearity (see, e.g., bicadze ).
Define the elliptic secondorder operator by the relation
(9) 
For (8), we have
By virtue of this, we obtain
Let satisfies the equation
(10) 
and the boundary conditions
(11) 
Under these conditions, for , we have the equation
(12) 
In view of (8), from (11), we obtain the following boundary condition:
(13) 
The equation (10) can be treated as a regularization of the Hamilton–Jacobi equation via the method of vanishing viscosity bookEvans . The boundary value problem (10), (11) produces an approximate solution of the problem (5), (6) for small values of :
In this case, is defined according to (8) from the solution of the linear boundary value problem (12), (13).
3 Numerical implementation
An approximate solution of the BVP (5)–(7) is represented (see (8)) in the form
(14) 
at a sufficiently low value of . In this case, is defined as the solution of the BVP (12), (13). In the present work, the numerical implementation of this approach is carried out on the basis of standard finiteelement approximations brenner2008mathematical ; larson2013finite . The main features of the computational algorithm result from the fact that the BVP of diffusion–reaction (12), (13) at small is singularly perturbed, i.e., we have a small parameter at higher derivatives holmes2012introduction ; verhulst2005methods .
Let us consider a standard quasiuniform triangulation of the domain into triangles in the 2D case or tetrahedra for 3D case. Let
Denote by and the linear finiteelement spaces.
For the BVP (12), (13), we put into the correspondence the variational problem of finding the numerical solution from the conditions
(15) 
By (9), for the bilinear form, we have
The differential problem (12), (13) satisfies the maximum principle. In particular, this guarantees the positiveness of the solution. More precisely (see, e.g., protter2012maximum ; gilbarg2015elliptic ), for points inside the domain , we have
This the most important property must be also fulfilled for the solution of the discrete problem (15):
(16) 
If (16) holds, we speak of monotone approximations for the solution of the diffusion–reaction problem.
Even for regular boundary value problems, where the parameter in (12) is not small, monotone approximations can be constructed using linear finite elements with restrictions on the computational grid (Delaunaytype mesh, see, for instance, letniowski1992three ; huang2011discrete ). Additional restrictions appear (see, e.g., ciarlet1973maximum ; brandts2008discrete ) on the magnitude of the reaction coefficient. With respect to our problem (12), (13), for the grid step size, we have .
Restrictions on the grid due to the reaction coefficient can be removed. The standard approach is related to the correction of approximations for the reaction coefficient based on the lumping procedure (see, e.g., Thomee2006 ).
The standard approach to the solution of singularly perturbed diffusion–reaction problems (see roos2008robust ; miller2012fitted ) is based on using computational grids with refinements in the vicinity of boundaries. A refinement of the grid is directly related to the value of the small parameter .
Another possibility to monotonize the solution of the problem (12), (13) at small values of is the following approach. As noted in the paper cai2014natural , for singularly perturbed problems for the diffusion–convection equation, the use of finiteelement approximations of higher order not only increases the accuracy of the approximate solution, but improves the monotonicity property as well. It is interesting to check whether there is the same effect in the numerical solution of singularly perturbed problems for the diffusion–reaction equations.
4 Numerical experiments
The 2D BVP (5)–(7) in the Lshaped region depicted in Fig. 1 is considered as a model problem. We start with the simplest case, when . The calculations have been performed on various grids. The basic (medium) computational grid, which contains 10,465 nodes and 20,480 triangles, is shown in Fig. 2.
In solving this problem, the key point is the dependence of the solution on the small parameter . The numerical solution obtained on a very fine grid with is treated as the exact one. The solution of the auxiliary problem (12), (13) under these conditions is presented in Fig. 3, and the corresponding function , determined according to (14), is shown in Fig. 4. The influence of the parameter can be observed in Fig. 5, where the solution in the cross section is plotted (the red line in Fig. 5). In our model problem, a good accuracy is achieved for .
The increase in accuracy can be achieved, first of all, by using finer grids. The solution for various on the coarse grid (2,673 nodes and 5,120 triangles) is given in Fig. 6. In this case, for , the solution is nonmonotone, i.e., at some nodes of the computational grid we have . Similar data for the basic grid are presented in Fig. 7. Here, the nonmonotonicity appears at . Figure 8 demonstrates the numerical results obtained on the fine grid (41,409 nodes and 81,920 triangles). The nonmonotonicity of the approximate solution occurs at .
In the practical use of the approach (12)–(14), it seems reasonable to follow the next strategy. We solve a number of auxiliary problems (12), (13) with a stepbystep decrease of the parameter as long as the maximum principle holds. The solution obtained with the smallest is taken as the approximate solution of the problem (5)–(7).
Our computational grids consist of rectangular isosceles triangles. Because of this, the nonmonotonicity is due to the reaction coefficient only. To monotonize discrete solutions, it is sufficient to apply the standard procedure of the reaction coefficient lumping Thomee2006 . The effect of diagonalization of the reactive term in the finiteelement approximation in predictions on different computational grids can be observed in Figures 9–11. In this case, the maximum principle holds for all .
The accuracy of the approximate solution decreases from some value of as the parameter decreases. Moreover, the value of this optimal value is close to the value that we had without the lumping procedure. Therefore, we can use the diagonalization procedure for selecting the parameter using the monotonicity condition for the discrete solution of the standard finiteelement approximation. In our case (see Figures 6–8), we select for the coarse grid, — for the basic grid and — for the fine grid.
Above, we have used linear finite elements. Below, we will present numerical results obtained on the basic grid for Lagrangian finite elements of degree , i.e., for approximations . Figure 12 demonstrates the approximate solution obtained using finite elements of degree 3. A comparison with the case (see Fig. 7) indicates that the solution is more accurate and, in addition, it is possible to carry out calculations with a smaller value of the parameter . These effects become more pronounced when using finite elements of degree 5 (see Fig. 13) and degree 7 (see Fig. 14).
Thus, the computational algorithm for solving the eikonal equation (BVP (5)–(7)) can be based on the solution of the auxiliary problem (12), (13). In doing so, we employ the minimum value of the parameter that provides the monotone solution on sufficiently fine computational grids using higher degree Lagrangian finite elements.
Special attention should be paid to the problem (5)–(7) in the anisotropic case. We have considered a variant with constant coefficients, where in (6), we had
The convergence of the approximate solution with decreasing is given in Fig. 15. Calculations have been performed on the basic grid using finiteelement approximation with . The numerical solution of the problem (5)–(7) for is shown in Fig. 16. Similar data obtained for a more pronounced anisotropy:
5 Conclusions

A Dirichlet problem is considered for the multidimensional eikonal equation in a bounded domain with an anisotropic medium. The main peculiarities of such problems results from the fact that the eikonal equation is nonlinear.

An approximate solution is constructed using a transformation of the original nonlinear boundary value problem to a linear boundary value problem for the diffusion–reaction equation for an auxiliary function. The transformed equation belongs to the class of singularly perturbed problems, i.e., there is a small parameter at higher derivatives.

Computational algorithms are constructed using standard finiteelement approximations on triangular (2D problems) or tetrahedral (3D problems) grids. Monotonization of a discrete solution is achieved not only by using finer grids, but also via a correction of approximations for the reaction coefficient using the lumping procedure. The use of finite elements of high degree is studied, too.

Numerical experiments have been performed for 2D problems in order to demonstrate the robustness of the approach proposed in the work for solving boundary value problems for the eikonal equation in an anisotropic medium. In particular, a good accuracy is observed when using sufficiently fine grids and Lagrangian finite elements of higher degree.
Acknowledgments
Petr Vabishchevich gratefully acknowledges support from the the Russian Federation Government (# 14.Y26.31.0013).
References
 (1) M. Born, E. Wolf, Principles of Optics: Electromagnetic Theory of Propagation, Interference and Diffraction of Light, Cambridge university press, 2005.
 (2) Y. A. Kravtsov, Y. I. Orlov, Geometrical Optics of Inhomogeneous Media, Springer, 1990.

(3)
J. A. Sethian, Level Set Methods and Fast Marching Methods: Evolving Interfaces in Computational Geometry, Fluid Mechanics, Computer Vision, and Materials Science, Vol. 3, Cambridge university press, 1999.
 (4) A. Gilles, K. Pierre, Mathematical Problems in Image Processing: Partial Differential Equations and the Calculus of Variations, Springer, 2006.
 (5) S. N. Kružkov, Generalized solutions of the HamiltonJacobi equations of eikonal type. I. Formulation of the problems; existence, uniqueness and stability theorems; some properties of the solutions, Sbornik: Mathematics 27 (3) (1975) 406–446.
 (6) P.L. Lions, Generalized Solutions of HamiltonJacobi Equations, Pitman, 1982.
 (7) J. Tsitsiklis, Fast marching methods, IEEE Transactions on Automatic Control 40 (1995) 1528–1538.
 (8) J. A. Sethian, Fast marching methods, SIAM Review 41 (2) (1999) 199–235.
 (9) Y.H. R. Tsai, L.T. Cheng, S. Osher, H.K. Zhao, Fast sweeping algorithms for a class of Hamilton–Jacobi equations, SIAM Journal on Numerical Analysis 41 (2) (2003) 673–694.
 (10) H. Zhao, Fast sweeping method for eikonal equations, Mathematics of Computation (74) (2005) 603–627.
 (11) W.K. Jeong, R. T. Whitaker, A fast iterative method for eikonal equations, SIAM Journal on Scientific Computing 30 (5) (2008) 2512–2534.
 (12) Z. Fu, W.K. Jeong, Y. Pan, R. M. Kirby, R. T. Whitaker, A fast iterative method for solving the eikonal equation on triangulated surfaces, SIAM Journal on Scientific Computing 33 (5) (2011) 2468–2488.
 (13) Z. Fu, R. M. Kirby, R. T. Whitaker, A fast iterative method for solving the eikonal equation on tetrahedral domains, SIAM Journal on Scientific Computing 35 (5) (2013) C473–C494.
 (14) J. V. Gómez, D. Alvarez, S. Garrido, L. Moreno, Fast methods for eikonal equations: an experimental survey, arXiv preprint arXiv:1506.03771.
 (15) A. G. Belyaev, P.A. Fayolle, On variational and PDEbased distance function approximations, in: Computer Graphics Forum, Vol. 34, Wiley Online Library, 2015, pp. 104–118.

(16)
C. Li, C. Xu, C. Gui, M. D. Fox, Level set evolution without reinitialization: a new variational formulation, in: Computer Vision and Pattern Recognition, 2005. CVPR 2005. IEEE Computer Society Conference on, Vol. 1, IEEE, 2005, pp. 430–436.
 (17) E. Fares, W. Schröder, A differential equation for approximate wall distance, International Journal for Numerical Methods in Fluids 39 (8) (2002) 743–762.
 (18) T. Bhattacharya, E. DiBenedetto, J. Manfredi, Limits as of and related extremal problems, Rendiconti del Seminario Matematico Universit‘a e Polytecnico di Torino 47 (1989) 15–68.
 (19) B. Kawohl, On a family of torsional creep problems, J. Reine Angew. Math 410 (1) (1990) 1–22.
 (20) K. S. Gurumoorthy, A. Rangarajan, A Schrödinger equation for the fast computation of approximate Euclidean distance functions, in: International Conference on Scale Space and Variational Methods in Computer Vision, Springer, 2009, pp. 100–111.
 (21) M. Sethi, A. Rangarajan, K. Gurumoorthy, The Schrödinger distance transform (SDT) for pointsets and curves, in: Computer Vision and Pattern Recognition (CVPR), 2012, IEEE, 2012, pp. 198–205.
 (22) K. Crane, C. Weischedel, M. Wardetzky, Geodesics in heat: A new approach to computing distance based on heat flow, ACM Transactions on Graphics (TOG) 32 (5) (2013) 152.
 (23) H. Roos, M. Stynes, L. Tobiska, Robust Numerical Methods for Singularly Perturbed Differential Equations: ConvectionDiffusionReaction and Flow Problems, Vol. 24, Springer Science & Business Media, 2008.

(24)
J. J. H. Miller, E. O’Riordan, G. I. Shishkin, Fitted Numerical Methods for Singular Perturbation Problems: Error Estimates in the Maximum Norm for Linear Problems in One and Two Dimensions, World Scientific, 2012.
 (25) S. C. Brenner, L. R. Scott, The Mathematical Theory of Finite Element Methods, Springer, New York, 2008.
 (26) M. G. Larson, F. Bengzon, The Finite Element Method: Theory, Implementation, and Applications, Springe, 2013.
 (27) A. V. Bitzadze, Some Classes of Partial Differential Equations, Nauka, 1981, in Russian.
 (28) L. C. Evans, Partial Differential Equations, AMS, 2010.
 (29) M. H. Holmes, Introduction to Perturbation Methods, Springer, 2012.
 (30) F. Verhulst, Methods and Applications of Singular Perturbations: Boundary Layers and Multiple Timescale Dynamics, Springe, 2005.
 (31) M. H. Protter, H. F. Weinberger, Maximum Principles in Differential Equations, Springer, 2012.
 (32) D. Gilbarg, N. S. Trudinger, Elliptic Partial Differential Equations of Second Order, springer, 2015.
 (33) F. W. Letniowski, Threedimensional Delaunay triangulations for finite element approximations to a secondorder diffusion operator, SIAM Journal on Scientific and Statistical Computing 13 (3) (1992) 765–770.
 (34) W. Huang, Discrete maximum principle and a delaunaytype mesh condition for linear finite element approximations of twodimensional anisotropic diffusion problems, Numerical Mathematics: Theory, Methods and Applications 4 (3) (2011) 319–334.
 (35) P.A. Ciarlet, P. G.and Raviart, Maximum principle and uniform convergence for the finite element method, Computer Methods in Applied Mechanics and Engineering 2 (1) (1973) 17–31.
 (36) J. H. Brandts, S. Korotov, M. Křížek, The discrete maximum principle for linear simplicial finite element approximations of a reaction–diffusion problem, Linear Algebra and its Applications 429 (10) (2008) 2344–2357.
 (37) V. Thomée, Galerkin Finite Element Methods for Parabolic Problems, Springer Verlag, Berlin, 2006.
 (38) Q. Cai, S. Kollmannsberger, E. SalaLardies, A. Huerta, E. Rank, On the natural stabilization of convection dominated problems using high order Bubnov–Galerkin finite elements, Computers & Mathematics with Applications 66 (12) (2014) 2545–2558.