Fast Marching with Higher Order Derivatives
This post is a short one and follows Escaping the maze: the Eikonal Equation unveiled, A numerical solution to the 1D eikonal equation and a Python implementation, The Godunov scheme for the solution of the two-dimensional eikonal equation and the Fast Sweeping Method and Two-dimensional Fast Marching.
In the above posts, numerical solutions to the 1D and 2D eikonal equations were given using first order approximations to the first order derivatives. In this post, we consider the use of second order approximations to the first order derivatives to improve the accuracy of the Fast Marching Method (FMM).
The codes associated to this post are available at this link.
The idea is simple. The eikonal equation can be discretized as:
Instead of using first order evaluations of first forward and backward derivatives, namely (for derivatives along the x-axis),
and
(same for derivatives along the y-axis), we use second order evaluations of first forward and backward derivatives, namely,
and
(same for derivatives along the y-axis), if certain conditions are met, otherwise a first order evaluation of the first order derivatives is employed.
The forward second order derivative along the x-axis is used if
(same for the forward second order derivative along the y-axis).
The backward second order derivative along the x-axis is used if
Here, known means that the corresponding value is converged and the node is in the upwind region (see above quoted previous post).
We now sketch the reason why switching from second order derivatives down to first order.
As already underlined in the above quoted posts, for the FMM, causality requires that the arrival time at any grid point should only depend on the arrival times at points that the wavefront has already visited. This ensures a physically consistent propagation of the wavefront, where information flows outward from the source and not backward. Second-order schemes may inadvertently use information from points where the wavefront has not yet arrived, violating causality since they do not strictly enforce the use of already-computed (upwind) values.
The new value of the eikonal at each grid point is computed solving the quadratic equation
Let us consider how the first term of the discretization of the eikonal equation contributes to the coefficients a, b and c. If backward second order derivatives along the x-axis are used, then
and
Similar expressions we have if forward second order derivatives are employed, possibly along the y-axis.
If, on the other side, backward first order derivatives are exploited, then
and
For points at the borders of the computational domain, first order derivatives are employed, backward or forward depending on the border. For example, if m=0, a forward first order derivative is adopted. On the other side, if m=M-1, backward first order derivative is used.
The solution procedure of the quadratic equation is similar to that of the first order FMM.
