A fast sweeping method for Eikonal equations

By Hongkai Zhao

Abstract

In this paper a fast sweeping method for computing the numerical solution of Eikonal equations on a rectangular grid is presented. The method is an iterative method which uses upwind difference for discretization and uses Gauss-Seidel iterations with alternating sweeping ordering to solve the discretized system. The crucial idea is that each sweeping ordering follows a family of characteristics of the corresponding Eikonal equation in a certain direction simultaneously. The method has an optimal complexity of upper O left-parenthesis upper N right-parenthesis for upper N grid points and is extremely simple to implement in any number of dimensions. Monotonicity and stability properties of the fast sweeping algorithm are proven. Convergence and error estimates of the algorithm for computing the distance function is studied in detail. It is shown that 2 Superscript n Gauss-Seidel iterations is enough for the distance function in n dimensions. An estimation of the number of iterations for general Eikonal equations is also studied. Numerical examples are used to verify the analysis.

1. Introduction

The Eikonal equation

StartLayout 1st Row with Label left-parenthesis 1.1 right-parenthesis EndLabel StartAbsoluteValue nabla u left-parenthesis bold x right-parenthesis EndAbsoluteValue equals f left-parenthesis bold x right-parenthesis comma bold x element-of upper R Superscript n Baseline comma EndLayout

with boundary condition, u left-parenthesis bold x right-parenthesis equals phi left-parenthesis bold x right-parenthesis comma bold x element-of normal upper Gamma subset-of upper R Superscript n , has many applications in optimal control, computer vision, geometric optics, path planing, etc. This nonlinear boundary value problem is a first order hyperbolic partial differential equation (PDE). Information propagates along characteristics from the boundary. Due to the nonlinearity, characteristics may intersect like the formation of shocks in hyperbolic conservation law. The solution is still continuous at these intersections but may not be differentiable. The existence and uniqueness of the viscosity solution are shown in Reference4.

There are two key ingredients for any numerical algorithm for the Eikonal equation. The first key ingredient is the derivation of a consistent and accurate discretization scheme, i.e., a numerical Hamiltonian. The numerical scheme has to follow the causality of the partial differential equation and has to deal with non-differentiability at intersections of characteristics properly. Since the problem is nonlinear, a large system of nonlinear equations has to be solved after the discretization. Hence the second key ingredient is an efficient method for solving the large nonlinear system.

There are mainly two types of approaches for solving the Eikonal equation. One approach is to transform it to a time dependent problem. For example, if we have u left-parenthesis bold x right-parenthesis equals 0 , bold x element-of normal upper Gamma , then u left-parenthesis bold x right-parenthesis is the first arrival time at bold x for a wave front starting at normal upper Gamma with a normal velocity that is equal to StartFraction 1 Over f left-parenthesis bold x right-parenthesis EndFraction . This can be solved by the level set method. In the control framework, a semi-Lagrangian scheme is obtained for Hamilton-Jacobi equations by discretizing in time the dynamic programming principle Reference9Reference10. However, many time steps may be needed for the convergence of the solution in the entire domain due to finite speed of propagation and CFL condition for time stability. The other approach is to treat the problem as a stationary boundary value problem and to design an efficient numerical algorithm to solve the system of nonlinear equations after discretization. For example, the fast marching method Reference19Reference16Reference11 is of this type. In the fast marching method, the update of the solution follows the causality in a sequential way; i.e., the solution is updated one grid point by one grid point in the order that the solution is strictly increasing (decreasing). Hence an upwind difference scheme and a heapsort algorithm are needed. The complexity is of order upper O left-parenthesis upper N log upper N right-parenthesis for upper N grid points, where the log upper N factor comes from the heapsort algorithm.

Here we present and analyze an iterative algorithm, called the fast sweeping method, for computing the numerical solution for the Eikonal equation on a rectangular grid in any number of space dimensions. The fast sweeping method is motivated by the work in Reference2 and was first used in Reference21 for computing the distance function. The main idea of the fast sweeping method is to use nonlinear upwind difference and Gauss-Seidel iterations with alternating sweeping ordering. In contrast to the fast marching method, the fast sweeping method follows the causality along characteristics in a parallel way; i.e., all characteristics are divided into a finite number of groups according to their directions and each Gauss-Seidel iteration with a specific sweeping ordering covers a group of characteristics simultaneously. The fast sweeping method is extremely simple to implement. The algorithm is optimal in the sense that a finite number of iterations is needed. So the complexity of the algorithm is upper O left-parenthesis upper N right-parenthesis for a total of upper N grid points. The number of iterations is independent of grid size. The accuracy is the same as any other method which solves the same system of discretized equations. The fast sweeping method has been extended to more general Hamilton-Jacobi equations Reference18Reference12. Extensions to high order discretization will be studied in future reports.

The idea of alternating sweeping ordering was also used in Danielsson’s algorithm Reference6. The algorithm computes the distance mapping, i.e., the relative left-parenthesis x comma y right-parenthesis coordinate of a grid point to its closest point using an iterative procedure. Danielsson’s algorithm is based on a strict dimension by dimension discrete formulation which in general does not follow the real characteristics of the distance function in two and higher dimensions and hence results in low accuracy and twice as many iterations compared to the fast sweeping method we present here. Danielsson’s algorithm does not work for distance functions to more general data sets such as the distance to a curve or a surface. Neither does it extend to general Eikonal equations. Recently another discrete approach that uses the idea of fast sweeping method was proposed in Reference17. It can compute the distance function more accurately but does not apply to general Eikonal equations either. Other related methods include a dynamic programming approach and post sweeping idea in Reference15 and a group marching method in Reference13.

Here is the outline of the paper. In Section 2 we present the scheme and the motivation behind it. We then show a few monotonicity properties and a maximum change principle for the fast sweeping algorithm in Section 3. In Section 4 we prove the convergence and error estimates for the distance function. We discuss the fast sweeping method for general Eikonal equations in Section 5. In Section 6 we present numerical results to verify our analysis and we show a few applications.

2. The fast sweeping algorithm and the motivation

We present the fast sweeping method for computing the viscosity solution u left-parenthesis bold x right-parenthesis greater-than-or-equal-to 0 for the model problem

StartLayout 1st Row with Label left-parenthesis 2.1 right-parenthesis EndLabel StartLayout 1st Row 1st Column StartAbsoluteValue nabla u left-parenthesis bold x right-parenthesis EndAbsoluteValue equals f left-parenthesis bold x right-parenthesis comma 2nd Column bold x element-of upper R Superscript n Baseline comma 2nd Row 1st Column u left-parenthesis bold x right-parenthesis equals 0 comma 2nd Column bold x element-of normal upper Gamma subset-of upper R Superscript n Baseline comma EndLayout EndLayout

where f left-parenthesis bold x right-parenthesis greater-than 0 . For simplicity the algorithm is presented in two dimensions. The extension to higher dimensions is straightforward. We use bold x Subscript i comma j to denote a grid point in the computational domain normal upper Omega , h to denote the grid size and u Subscript i comma j Superscript h to denote the numerical solution at bold x Subscript i comma j .

2.1. The fast sweeping algorithm

Discretization

We use the following Godunov upwind difference scheme Reference14 to discretize the partial differential equation at interior grid points: StartLayout 1st Row with Label left-parenthesis 2.2 right-parenthesis EndLabel StartLayout 1st Row left-bracket left-parenthesis u Subscript i comma j Superscript h Baseline minus u Subscript x min Superscript h Baseline right-parenthesis Superscript plus Baseline right-bracket squared plus left-bracket left-parenthesis u Subscript i comma j Superscript h Baseline minus u Subscript y min Superscript h Baseline right-parenthesis Superscript plus Baseline right-bracket squared equals f Subscript i comma j Superscript 2 Baseline h squared 2nd Row Blank 3rd Row i equals 2 comma period period period comma upper I minus 1 comma j equals 2 comma period period period comma upper J minus 1 comma EndLayout EndLayout

where u Subscript x min Superscript h Baseline equals min left-parenthesis u Subscript i minus 1 comma j Superscript h Baseline comma u Subscript i plus 1 comma j Superscript h Baseline right-parenthesis comma u Subscript y min Superscript h Baseline equals min left-parenthesis u Subscript i comma j minus 1 Superscript h Baseline comma u Subscript i comma j plus 1 Superscript h Baseline right-parenthesis and left-parenthesis x right-parenthesis Superscript plus Baseline equals StartLayout Enlarged left-brace 1st Row 1st Column x comma 2nd Column x greater-than 0 comma 2nd Row 1st Column 0 comma 2nd Column x less-than-or-equal-to 0 period EndLayout

One sided difference is used at the boundary of the computational domain. For example, at a left boundary point bold x Subscript 1 comma j , a one sided difference is used in the x direction, left-bracket left-parenthesis u Subscript 1 comma j Superscript h Baseline minus u Subscript 2 comma j Superscript h Baseline right-parenthesis Superscript plus Baseline right-bracket squared plus left-bracket left-parenthesis u Subscript 1 comma j Superscript h Baseline minus u Subscript y min Superscript h Baseline right-parenthesis Superscript plus Baseline right-bracket squared equals f Subscript 1 comma j Superscript 2 Baseline h squared period

Initialization

To enforce the boundary condition, u left-parenthesis bold x right-parenthesis equals 0 for bold x element-of normal upper Gamma subset-of upper R Superscript n , assign exact values or interpolated values at grid points in or near normal upper Gamma . These values are fixed in later calculations. Assign large positive values at all other grid points. These values will be updated later.

Gauss-Seidel iterations with alternating sweeping orderings

At each grid bold x Subscript i comma j whose value is not fixed during the initialization, compute the solution, denoted by u overbar , of (Equation2.2) from the current values of its neighbors u Subscript i plus-or-minus 1 comma j Superscript h Baseline comma u Subscript i comma j plus-or-minus 1 Superscript h and then update u Subscript i comma j Superscript h to be the smaller one between u overbar and its current value; i.e., u Subscript i comma j Superscript n e w Baseline equals min left-parenthesis u Subscript i comma j Superscript o l d Baseline comma u overbar right-parenthesis . We sweep the whole domain with four alternating orderings repeatedly: StartLayout 1st Row 1st Column left-parenthesis 1 right-parenthesis i equals 1 colon upper I comma j equals 1 colon upper J comma 2nd Column left-parenthesis 2 right-parenthesis i equals upper I colon 1 comma j equals 1 colon upper J comma 2nd Row 1st Column left-parenthesis 3 right-parenthesis i equals upper I colon 1 comma j equals upper J colon 1 comma 2nd Column left-parenthesis 4 right-parenthesis i equals 1 colon upper I comma j equals upper J colon 1 period EndLayout

The unique solution to the equation

StartLayout 1st Row with Label left-parenthesis 2.3 right-parenthesis EndLabel left-bracket left-parenthesis x minus a right-parenthesis Superscript plus Baseline right-bracket squared plus left-bracket left-parenthesis x minus b right-parenthesis Superscript plus Baseline right-bracket squared equals f Subscript i comma j Superscript 2 Baseline h squared comma EndLayout

where a equals u Subscript x min Superscript h Baseline comma b equals u Subscript y min Superscript h , is

StartLayout 1st Row with Label left-parenthesis 2.4 right-parenthesis EndLabel x overbar equals StartLayout Enlarged left-brace 1st Row 1st Column min left-parenthesis a comma b right-parenthesis plus f Subscript i comma j Baseline h comma 2nd Column StartAbsoluteValue a minus b EndAbsoluteValue greater-than-or-equal-to f Subscript i comma j Baseline h comma 2nd Row 1st Column StartFraction a plus b plus StartRoot 2 f Subscript i comma j Superscript 2 Baseline h squared minus left-parenthesis a minus b right-parenthesis squared EndRoot Over 2 EndFraction comma 2nd Column StartAbsoluteValue a minus b EndAbsoluteValue less-than f Subscript i comma j Baseline h period EndLayout EndLayout

In n dimensions the unique solution x overbar to

StartLayout 1st Row with Label left-parenthesis 2.5 right-parenthesis EndLabel left-bracket left-parenthesis x minus a 1 right-parenthesis Superscript plus Baseline right-bracket squared plus left-bracket left-parenthesis x minus a 2 right-parenthesis Superscript plus Baseline right-bracket squared plus midline-horizontal-ellipsis plus left-bracket left-parenthesis x minus a Subscript n Baseline right-parenthesis Superscript plus Baseline right-bracket squared equals f Subscript i comma j Superscript 2 Baseline h squared EndLayout

can be found in the following systematic way. First we order the a Subscript k ’s in increasing order. Without loss of generality, assume a 1 less-than-or-equal-to a 2 less-than-or-equal-to midline-horizontal-ellipsis less-than-or-equal-to a Subscript n and define a Subscript n plus 1 Baseline equals normal infinity . There is an integer p , 1 less-than-or-equal-to p less-than-or-equal-to n , such that x overbar is the unique solution that satisfies

StartLayout 1st Row with Label left-parenthesis 2.6 right-parenthesis EndLabel left-parenthesis x minus a 1 right-parenthesis squared plus left-parenthesis x minus a 2 right-parenthesis squared plus midline-horizontal-ellipsis plus left-parenthesis x minus a Subscript p Baseline right-parenthesis squared equals f Subscript i comma j Superscript 2 Baseline h squared and a Subscript p Baseline less-than x overbar less-than-or-equal-to a Subscript p plus 1 Baseline semicolon EndLayout

i.e., x overbar is the intersection of the straight line bold x equals bold y comma bold x comma bold y element-of upper R Superscript p , with the sphere centered at bold a equals left-parenthesis a 1 comma a 2 comma period period period comma a Subscript p Baseline right-parenthesis of radius f Subscript i comma j Baseline h in the first quadrant in upper R Superscript p . We find x overbar and p in the following recursive way. Start with p equals 1 . If x overTilde equals a 1 plus f Subscript i comma j Baseline h less-than-or-equal-to a 2 , then x overbar equals x overTilde . Otherwise find the unique solution x overTilde greater-than a 2 that satisfies

left-parenthesis x minus a 1 right-parenthesis squared plus left-parenthesis x minus a 2 right-parenthesis squared equals f Subscript i comma j Superscript 2 Baseline h squared period

If x overTilde less-than-or-equal-to a 3 , then x overbar equals x overTilde . Otherwise repeat the procedure until we find p and x overbar that satisfy (Equation2.6).

Here are a few remarks about the algorithm:

(1)

The upwind difference scheme (Equation2.2) is a special case of the general Godunov numerical Hamiltonian proposed in Reference1. The numerical Hamiltonian can be written as StartLayout 1st Row 1st Column Blank 2nd Column upper H Superscript h Baseline left-parenthesis upper D Subscript minus Superscript x Baseline u Subscript i comma j Baseline comma upper D Subscript plus Superscript x Baseline u Subscript i comma j Baseline comma upper D Subscript minus Superscript y Baseline u Subscript i comma j Baseline comma upper D Subscript plus Superscript y Baseline u Subscript i comma j Baseline right-parenthesis 2nd Row 1st Column Blank 2nd Column equals StartRoot max StartSet left-parenthesis upper D Subscript minus Superscript x Baseline u Subscript i comma j Baseline right-parenthesis Superscript plus Baseline comma left-parenthesis upper D Subscript plus Superscript x Baseline u Subscript i comma j Baseline right-parenthesis Superscript minus Baseline EndSet squared plus max StartSet left-parenthesis upper D Subscript minus Superscript y Baseline u Subscript i comma j Baseline right-parenthesis Superscript plus Baseline comma left-parenthesis upper D Subscript plus Superscript y Baseline u Subscript i comma j Baseline right-parenthesis Superscript minus Baseline EndSet squared EndRoot comma EndLayout

where StartLayout 1st Row upper D Subscript minus Superscript x Baseline u Subscript i comma j Baseline equals StartFraction u Subscript i comma j Baseline minus u Subscript i minus 1 comma j Baseline Over h EndFraction comma upper D Subscript plus Superscript x Baseline u Subscript i comma j Baseline equals StartFraction u Subscript i plus 1 comma j Baseline minus u Subscript i comma j Baseline Over h EndFraction comma 2nd Row upper D Subscript minus Superscript y Baseline u Subscript i comma j Baseline equals StartFraction u Subscript i comma j Baseline minus u Subscript i comma j minus 1 Baseline Over h EndFraction comma upper D Subscript plus Superscript y Baseline u Subscript i comma j Baseline equals StartFraction u Subscript i comma j plus 1 Baseline minus u Subscript i comma j Baseline Over h EndFraction comma EndLayout

and left-parenthesis dot right-parenthesis Superscript plus means taking the positive part and left-parenthesis dot right-parenthesis Superscript minus means taking the negative part. Instead of using the upwind difference scheme (Equation2.2), we can also use StartLayout 1st Row with Label left-parenthesis 2.7 right-parenthesis EndLabel StartLayout 1st Row 1st Column left-bracket left-parenthesis u Subscript i comma j Superscript h Baseline 2nd Column minus u Subscript i minus 1 comma j Superscript h Baseline right-parenthesis Superscript plus Baseline right-bracket Superscript 2 Baseline plus left-bracket left-parenthesis u Subscript i comma j Superscript h Baseline minus u Subscript i plus 1 comma j Superscript h Baseline right-parenthesis Superscript plus Baseline right-bracket squared 2nd Row 1st Column Blank 2nd Column plus left-bracket left-parenthesis u Subscript i comma j Superscript h Baseline minus u Subscript i comma j minus 1 Superscript h Baseline right-parenthesis Superscript plus Baseline right-bracket squared plus left-bracket left-parenthesis u Subscript i comma j Superscript h Baseline minus u Subscript i comma j plus 1 Superscript h Baseline right-parenthesis Superscript plus Baseline right-bracket squared equals f Subscript i comma j Superscript 2 Baseline h squared comma 3rd Row 1st Column Blank 2nd Column i equals 2 comma period period period comma upper I minus 1 comma j equals 2 comma period period period comma upper J minus 1 comma EndLayout EndLayout

which corresponds to the numerical Hamiltonian StartLayout 1st Row 1st Column Blank 2nd Column upper H Superscript h Baseline left-parenthesis upper D Subscript minus Superscript x Baseline u Subscript i comma j Baseline comma upper D Subscript plus Superscript x Baseline u Subscript i comma j Baseline comma upper D Subscript minus Superscript y Baseline u Subscript i comma j Baseline comma upper D Subscript plus Superscript y Baseline u Subscript i comma j Baseline right-parenthesis 2nd Row 1st Column Blank 2nd Column equals StartRoot left-bracket left-parenthesis upper D Subscript minus Superscript x Baseline u Subscript i comma j Baseline right-parenthesis Superscript plus Baseline right-bracket squared plus left-bracket left-parenthesis upper D Subscript plus Superscript x Baseline u Subscript i comma j Baseline right-parenthesis Superscript minus Baseline right-bracket squared plus left-bracket left-parenthesis upper D Subscript minus Superscript y Baseline u Subscript i comma j Baseline right-parenthesis Superscript plus Baseline right-bracket squared plus left-bracket left-parenthesis upper D Subscript plus Superscript y Baseline u Subscript i comma j Baseline right-parenthesis Superscript minus Baseline right-bracket squared EndRoot period EndLayout

Both numerical Hamiltonians are monotone. The only difference between these two formulations is at the intersections of characteristics. The second formulation may have larger truncation errors at those intersections as will be explained in Section 4.

(2)

In practice it may be desirable to restrict the computation to a neighborhood of the boundary normal upper Gamma . For example, if we want to restrict the computation in the neighborhood where the first arrival time is less than upper T , i.e., StartSet bold x Subscript i comma j Baseline colon u left-parenthesis bold x Subscript i comma j Baseline right-parenthesis less-than upper T EndSet , then we can use the following simple cutoff criterion: in the Gauss-Seidel iteration we update the solution at a grid point bold x Subscript i comma j only if at least one of its neighbors has a value smaller than upper T , i.e., if min left-parenthesis u Subscript x min Superscript h Baseline comma u Subscript y min Superscript h Baseline right-parenthesis less-than upper T .

(3)

The large value assigned initially should be larger than the maximum possible value of u left-parenthesis bold x right-parenthesis in the computation domain. For example, let f Subscript upper M Baseline equals max Underscript bold x element-of normal upper Omega Endscripts f left-parenthesis bold x right-parenthesis and let upper D be the diameter of the computational domain normal upper Omega . The initially assigned large value should be larger than upper F Subscript upper M Baseline upper D .

(4)

In higher dimensions, the discretization (Equation2.2) is easily extended dimension by dimension and there are 2 Superscript n different sweeping orderings in n dimensions.

(5)

If we want to compute the viscosity solution u left-parenthesis bold x right-parenthesis less-than-or-equal-to 0 for (Equation2.1), we modify the discretization (Equation2.2) to StartLayout 1st Row with Label left-parenthesis 2.8 right-parenthesis EndLabel StartLayout 1st Row left-bracket left-parenthesis u Subscript i comma j Superscript h Baseline minus u Subscript x max Superscript h Baseline right-parenthesis Superscript minus Baseline right-bracket squared plus left-bracket left-parenthesis u Subscript i comma j Superscript h Baseline minus u Subscript y max Superscript h Baseline right-parenthesis Superscript minus Baseline right-bracket squared equals f Subscript i comma j Superscript 2 Baseline h squared comma 2nd Row i equals 2 comma period period period comma upper I minus 1 comma j equals 2 comma period period period comma upper J minus 1 comma EndLayout EndLayout

where u Subscript x max Superscript h Baseline equals max left-parenthesis u Subscript i minus 1 comma j Superscript h Baseline comma u Subscript i plus 1 comma j Superscript h Baseline right-parenthesis comma u Subscript y max Superscript h Baseline equals max left-parenthesis u Subscript i comma j minus 1 Superscript h Baseline comma u Subscript i comma j plus 1 Superscript h Baseline right-parenthesis and left-parenthesis x right-parenthesis Superscript minus Baseline equals StartLayout Enlarged left-brace 1st Row 1st Column 0 comma 2nd Column x greater-than 0 comma 2nd Row 1st Column negative x comma 2nd Column x less-than-or-equal-to 0 period EndLayout

In the initialization, we assign small negative values at grid points whose values are to be updated. In the Gauss-Seidel iteration, we update the value at a grid point only if the new value obtained by solving (Equation2.8) is larger than its old value.

2.2. The motivation

In the fast sweeping algorithm the upwind difference scheme used in the discretization enforces the causality; i.e., the solution at a grid point is determined by its neighboring values that are smaller. The one sided difference scheme at the boundary enforces the propagation of information to be from inside to outside since the data set normal upper Gamma is contained in the computational domain. If all grid points can be ordered according to the causality along characteristics, one iteration of the Gauss-Seidel iteration is enough for convergence. For example, the heapsort algorithm is used in the fast marching method to sort out this order every time a grid point is updated. The key point behind Gauss-Seidel iterations with different sweeping ordering is that each sweep will follow the causality of a group of characteristics in certain directions simultaneously and all characteristics can be divided into a finite number of such groups according to their directions. The value at each grid point is always nonincreasing during the iterations due to the updating rule. Whenever a grid point obtains the minimal value it can reach, the value is the correct value and the value will not be changed in later iterations.

We use the distance function as an example to illustrate the motivation. The distance function d left-parenthesis bold x right-parenthesis to a set normal upper Gamma satisfies the Eikonal equation

StartAbsoluteValue nabla d left-parenthesis bold x right-parenthesis EndAbsoluteValue equals 1 comma d left-parenthesis bold x right-parenthesis equals 0 comma bold x element-of normal upper Gamma period

All characteristics of this equation are straight lines that radiate from the set normal upper Gamma . In one dimension, the upwind differencing at the interior grid point i is

StartLayout 1st Row with Label left-parenthesis 2.9 right-parenthesis EndLabel left-bracket left-parenthesis u Subscript i Superscript h Baseline minus min left-parenthesis u Subscript i minus 1 Superscript h Baseline comma u Subscript i plus 1 Superscript h Baseline right-parenthesis right-parenthesis Superscript plus Baseline right-bracket squared equals h squared comma 2 less-than-or-equal-to i less-than-or-equal-to upper I minus 1 period EndLayout

We use two Gauss-Seidel iterations with sweeping orderings, i equals 1 colon upper I and i equals upper I colon 1 successively, to solve the above system. The update of the distance value at grid i simply becomes

u Subscript i Superscript n e w Baseline equals min left-parenthesis min left-parenthesis u Subscript i minus 1 Baseline comma u Subscript i plus 1 Baseline right-parenthesis plus h comma u Subscript i Baseline right-parenthesis period

Figure 2.1 shows how one sweep from left to right followed by one more sweep from right to left is enough to finish the calculation of the distance function. This follows because there are only two directions for the characteristics in one dimension, left to right or vice versa. In another word, the distance value at any grid point can be computed from either its left neighbor or right neighbor by exactly d Subscript i Baseline equals min left-parenthesis d Subscript i minus 1 Baseline comma d Subscript i plus 1 Baseline right-parenthesis plus h . The first sweep will cover those characteristics that go from left to right; i.e., those grid points whose values are determined by their left neighbors are computed correctly. Similarly, in the second sweep all those grid points whose values are determined by their right neighbors are computed correctly. Since we only update the current value if the newly computed value is smaller, those values that have been calculated correctly in the first sweep have achieved their minimal possible values and will not be changed in the second sweep. Convergence in two sweeps is true for arbitrary Eikonal equations in one dimension. In the special case of the distance function, it is easy to see that the fast sweeping method finds the exact distance function in two sweeps.

In higher dimensions characteristics have infinitely many directions which cannot be followed exactly by the Cartesian grid lines. Here are two important questions for the fast sweeping algorithm: (1) How many Gauss-Seidel iterations are needed? (2) What is the error estimate? The most important observation is that all directions of characteristics can be classified into a finite number of groups for distance functions. For example, in two dimensions all directions of characteristics can be classified into four groups, up-right, up-left, down-left and down-right. Information propagates along characteristics in the above four groups of directions. The four different orderings of the Gauss-Seidel iterations and the upwind differencing are meant to cover the four groups of characteristics, respectively. Figure 2.2(a) illustrates why the fast sweeping method converges after four sweeps with different orderings for the distance function to a single point. The solution u Subscript i comma j Superscript h at each grid point in the first quadrant depends on the solution at u Subscript i minus 1 comma j Superscript h and u Subscript i comma j minus 1 Superscript h which have already been computed and can be recursively traced all the way back to the data point in the first sweep. So we get the correct values for all grid points in the first quadrant plus the points on the positive x and y axes after the first sweep. For the same reason, after the second sweep the grid points in the second quadrant and on the negative x axis get the correct values. Moreover, since those grid points in the first quadrant and on the positive x and y axes already have their minimal values and satisfy the discretized equations, these values will not change in the second sweep. Similarly, after the third sweep those grid points in the third quadrant and on the negative y axis get the correct values and the computed correct values in the first and second quadrants are maintained. After four sweeps we get the correct values for all grid points that satisfy the system of equations (Equation2.2). Figure 2.2(b) demonstrates another case which computes the distance function to a circle in two dimensions using the fast sweeping algorithm. Again, each grid point gets its correct value in one of the four sweeps.

In the case of one data point, the distance function is smooth except at the data point. For a more general data set, interactions of characteristics at their intersections can cause more than 2 Superscript n sweeps for the iteration to converge in n dimensions. It is impossible to track the exact number of sweeps for the highly nonlinear discretized system in general. However, we will show that in n dimensions, after 2 Superscript n sweeps the fast sweeping method can compute a numerical solution to the distance function that is as accurate as the numerical solution after the iteration converges. This means 2 Superscript n sweeps is good enough in practice for computing the distance function to an arbitrary data set. For general Eikonal equations, the characteristics are curves instead of straight lines. So more than one sweep may be needed to cover one characteristic curve. We will see that given a fixed domain and the right-hand side f left-parenthesis bold x right-parenthesis , the number of sweeps needed is still finite and is independent of grid size.

3. Basic properties of the fast sweeping algorithm

Here we prove a few basic monotonicity properties and a maximum change principle for the fast sweeping algorithm. The following simple fact for the solution of equation (Equation2.5) provides the monotonicity and maximum change principle for the fast sweeping algorithm.

Lemma 3.1

Let x overbar be the solution to equation Equation2.5. We have

StartLayout 1st Row with Label left-parenthesis 3.1 right-parenthesis EndLabel 1 greater-than-or-equal-to StartFraction partial-differential x overbar Over partial-differential a 1 EndFraction greater-than-or-equal-to StartFraction partial-differential x overbar Over partial-differential a 2 EndFraction greater-than-or-equal-to midline-horizontal-ellipsis greater-than-or-equal-to StartFraction partial-differential x overbar Over partial-differential a Subscript n Baseline EndFraction greater-than-or-equal-to 0 EndLayout

and

StartLayout 1st Row with Label left-parenthesis 3.2 right-parenthesis EndLabel StartFraction partial-differential x overbar Over partial-differential a 1 EndFraction plus StartFraction partial-differential x overbar Over partial-differential a 2 EndFraction plus midline-horizontal-ellipsis plus StartFraction partial-differential x overbar Over partial-differential a Subscript n Baseline EndFraction equals 1 period EndLayout

Proof.

Differentiating equation (Equation2.5) with respect to a Subscript k , we get

StartFraction partial-differential x overbar Over partial-differential a Subscript k Baseline EndFraction equals StartFraction left-parenthesis x overbar minus a Subscript k Baseline right-parenthesis Superscript plus Baseline Over sigma-summation Underscript m equals 1 Overscript n Endscripts left-parenthesis x overbar minus a Subscript m Baseline right-parenthesis Superscript plus Baseline EndFraction period

As a direct consequence of the lemma, we have

Proposition 3.2

In the Gauss-Seidel iteration for the fast sweeping method, the maximum change of u Superscript h at any grid point is less than or equal to the maximum change of u Superscript h at its neighboring points.

Lemma 3.3

The fast sweeping algorithm is monotone in the initial data.

Proof.

This is a direct consequence from Lemma 3.1; i.e., the monotonicity property of the solution to (Equation2.5). If u Superscript h Baseline left-parenthesis bold x Subscript i comma j Baseline right-parenthesis less-than-or-equal-to v Superscript h Baseline left-parenthesis bold x Subscript i comma j Baseline right-parenthesis at all grid points initially, then u Superscript h Baseline left-parenthesis bold x Subscript i comma j Baseline right-parenthesis less-than-or-equal-to v Superscript h Baseline left-parenthesis bold x Subscript i comma j Baseline right-parenthesis at all grid points after any number of Gauss-Seidel iterations.

Lemma 3.4

The solution of the fast sweeping algorithm is nonincreasing with each Gauss-Seidel iteration.

Proof.

This is exactly because of the way we update the solution in the Gauss-Seidel iterations described in the third step of the algorithm.

The following corollary shows the stability property for the fast sweeping method. We present it in two dimensions but it is true in general.

Corollary 3.5

Let u Superscript left-parenthesis k right-parenthesis and v Superscript left-parenthesis k right-parenthesis be two numerical solutions at the k -th iteration of the fast sweeping algorithm. Let double-vertical-bar dot double-vertical-bar Subscript normal infinity be the maximum norm. We have

(1) double-vertical-bar u Superscript left-parenthesis k right-parenthesis Baseline minus v Superscript left-parenthesis k right-parenthesis Baseline double-vertical-bar Subscript normal infinity Baseline less-than-or-equal-to double-vertical-bar u Superscript left-parenthesis k minus 1 right-parenthesis Baseline minus v Superscript left-parenthesis k minus 1 right-parenthesis Baseline double-vertical-bar Subscript normal infinity ,

(2) 0 less-than-or-equal-to max Underscript i comma j Endscripts left-brace u Subscript i comma j Superscript left-parenthesis k right-parenthesis Baseline minus u Subscript i comma j Superscript left-parenthesis k plus 1 right-parenthesis Baseline right-brace less-than-or-equal-to max Underscript i comma j Endscripts left-brace u Subscript i comma j Superscript left-parenthesis k minus 1 right-parenthesis Baseline minus u Subscript i comma j Superscript left-parenthesis k right-parenthesis Baseline right-brace .

Proof.

Let us assume that the first update at the k -th iteration is at point bold x Subscript i comma j , u Subscript i comma j Superscript left-parenthesis k right-parenthesis Baseline equals min left-brace u Subscript i comma j Superscript left-parenthesis k minus 1 right-parenthesis Baseline comma u overbar right-brace , where u overbar solves (Equation2.2) with neighboring values u Subscript i minus 1 comma j Superscript left-parenthesis k minus 1 right-parenthesis , u Subscript i plus 1 comma j Superscript left-parenthesis k minus 1 right-parenthesis , u Subscript i comma j minus 1 Superscript left-parenthesis k minus 1 right-parenthesis , u Subscript i comma j plus 1 Superscript left-parenthesis k minus 1 right-parenthesis . The same is true for v Subscript i comma j Superscript left-parenthesis k right-parenthesis . From the maximum change principle, we have

StartAbsoluteValue u Subscript i comma j Superscript left-parenthesis k right-parenthesis Baseline minus v Subscript i comma j Superscript left-parenthesis k right-parenthesis Baseline EndAbsoluteValue less-than-or-equal-to double-vertical-bar u Superscript left-parenthesis k minus 1 right-parenthesis Baseline minus v Superscript left-parenthesis k minus 1 right-parenthesis Baseline double-vertical-bar Subscript normal infinity Baseline period

For an update at any other grid point later in the iteration, the neighboring values used for the update are either from the previous iteration or from an earlier update in the current iteration, both of which satisfy the above bound. By induction, we prove (1).

The second statement is a simple consequence of the monotonicity of the fast sweeping method and the previous statement by setting v Superscript left-parenthesis k right-parenthesis Baseline equals u Superscript left-parenthesis k minus 1 right-parenthesis .

Remark

The second statement provides an effective stopping criterion for sweeping. Before correct information from the boundary normal upper Gamma has reached all grid points, double-vertical-bar u Superscript left-parenthesis k right-parenthesis Baseline minus u Superscript left-parenthesis k minus 1 right-parenthesis Baseline double-vertical-bar Subscript normal infinity is of upper O left-parenthesis 1 right-parenthesis . When double-vertical-bar u Superscript left-parenthesis k right-parenthesis Baseline minus u Superscript left-parenthesis k minus 1 right-parenthesis Baseline double-vertical-bar Subscript normal infinity is upper O left-parenthesis h right-parenthesis , the information has reached all points and we can stop the sweeping in practice. Although the iteration has not converged yet, the numerical solution is already as accurate as the converged one as we will see from the numerical examples. The iterative solution will change less and less in later iterations and converges to the solution of the discretized system.

Theorem 3.6

The iterative solution by the fast sweeping algorithm converges monotonically to the solution of the discretized system.

Proof.

Denote the numerical solution after the k -th sweep by u Subscript i comma j Superscript left-parenthesis k right-parenthesis . Since u Subscript i comma j Superscript left-parenthesis k right-parenthesis is bounded below by 0 and is nonincreasing with Gauss-Seidel iterations, u Subscript i comma j Superscript left-parenthesis k right-parenthesis is convergent for all i comma j . After each sweep, for each i comma j we have

left-bracket left-parenthesis u Subscript i comma j Superscript left-parenthesis k right-parenthesis Baseline minus u Subscript x min Superscript left-parenthesis k right-parenthesis Baseline right-parenthesis Superscript plus Baseline right-bracket squared plus left-bracket left-parenthesis u Subscript i comma j Superscript left-parenthesis k right-parenthesis Baseline minus u Subscript y min Superscript left-parenthesis k right-parenthesis Baseline right-parenthesis Superscript plus Baseline right-bracket squared minus f Subscript i comma j Superscript 2 Baseline h squared greater-than-or-equal-to 0 comma

because any later update of neighbors of u Subscript i comma j Superscript left-parenthesis k right-parenthesis in the same sweep is nonincreasing. Moreover, it is easy to see that after u Subscript i comma j Superscript left-parenthesis k right-parenthesis is updated, the function

upper F left-parenthesis u Subscript i minus 1 comma j Superscript left-parenthesis k right-parenthesis Baseline comma u Subscript i plus 1 comma j Superscript left-parenthesis k right-parenthesis Baseline comma u Subscript i comma j minus 1 Superscript left-parenthesis k right-parenthesis Baseline comma u Subscript i comma j plus 1 Superscript left-parenthesis k right-parenthesis Baseline right-parenthesis equals left-bracket left-parenthesis u Subscript i comma j Superscript left-parenthesis k right-parenthesis Baseline minus u Subscript x min Superscript left-parenthesis k right-parenthesis Baseline right-parenthesis Superscript plus Baseline right-bracket squared plus left-bracket left-parenthesis u Subscript i comma j Superscript left-parenthesis k right-parenthesis Baseline minus u Subscript y min Superscript left-parenthesis k right-parenthesis Baseline right-parenthesis Superscript plus Baseline right-bracket squared minus f Subscript i comma j Superscript 2 Baseline h squared

is Lipschitz continuous in all variables and the Lipschitz constant is bounded by

2 max left-brace left-parenthesis u Subscript i comma j Superscript left-parenthesis k right-parenthesis Baseline minus u Subscript i minus 1 comma j Superscript left-parenthesis k right-parenthesis Baseline right-parenthesis Superscript plus Baseline comma left-parenthesis u Subscript i comma j Superscript left-parenthesis k right-parenthesis Baseline minus u Subscript i plus 1 comma j Superscript left-parenthesis k right-parenthesis Baseline right-parenthesis Superscript plus Baseline comma left-parenthesis u Subscript i comma j Superscript left-parenthesis k right-parenthesis Baseline minus u Subscript i comma j minus 1 Superscript left-parenthesis k right-parenthesis Baseline right-parenthesis Superscript plus Baseline comma left-parenthesis u Subscript i comma j Superscript left-parenthesis k right-parenthesis Baseline minus u Subscript i comma j plus 1 Superscript left-parenthesis k right-parenthesis Baseline right-parenthesis Superscript plus Baseline right-brace period

Since u Subscript i comma j Superscript left-parenthesis k right-parenthesis is monotonically convergent for every i comma j , we can have an upper bound upper C greater-than 0 for the Lipschitz constant. Let delta Superscript left-parenthesis k right-parenthesis Baseline equals max Underscript i comma j Endscripts left-parenthesis u Subscript i comma j Superscript left-parenthesis k minus 1 right-parenthesis Baseline minus u Subscript i comma j Superscript left-parenthesis k right-parenthesis Baseline right-parenthesis be the maximum change at all grid points during the k -th sweep. From Corollary 3.5 and the convergence of u Subscript i comma j Superscript left-parenthesis k right-parenthesis , delta Superscript left-parenthesis k right-parenthesis goes monotonically to zero. After the k -th iteration, we have

0 less-than-or-equal-to left-bracket left-parenthesis u Subscript i comma j Superscript left-parenthesis k right-parenthesis Baseline minus u Subscript x min Superscript left-parenthesis k right-parenthesis Baseline right-parenthesis Superscript plus Baseline right-bracket squared plus left-bracket left-parenthesis u Subscript i comma j Superscript left-parenthesis k right-parenthesis Baseline minus u Subscript y min Superscript left-parenthesis k right-parenthesis Baseline right-parenthesis Superscript plus Baseline right-bracket squared minus f Subscript i comma j Superscript 2 Baseline h squared less-than-or-equal-to upper C delta Superscript left-parenthesis k right-parenthesis Baseline period

So u Subscript i comma j Superscript left-parenthesis k right-parenthesis converges to the solution to (Equation2.2).

4. Convergence and error estimate for the distance function

Although it is shown that the fast sweeping algorithm is convergent by Theorem 3.6, the main issue for the efficiency of the fast sweeping method is the number of iterations needed and the error estimate. In this section we prove a few concrete results for the distance function. Since the fast sweeping algorithm is highly nonlinear, the proof can only be based on the monotonicity properties and the maximum change principle proved in Section 3. Some of the proofs will be stated in two dimensions for simplicity. We use the following notations: normal upper Gamma denotes the data set to which we want to compute the distance function, u Superscript h Baseline left-parenthesis bold x comma normal upper Gamma right-parenthesis denotes the numerical solution using the fast sweeping method, d left-parenthesis bold x comma normal upper Gamma right-parenthesis denotes the exact distance function, h is the grid size, and n is the spatial dimension.

Theorem 4.1

For a single data point normal upper Gamma equals StartSet bold x 0 EndSet , the numerical solution, u Superscript h Baseline left-parenthesis bold x comma bold x 0 right-parenthesis , of the fast sweeping method converges in 2 Superscript n sweeps in upper R Superscript n and satisfies

d left-parenthesis bold x comma bold x 0 right-parenthesis less-than-or-equal-to u Superscript h Baseline left-parenthesis bold x comma bold x 0 right-parenthesis less-than-or-equal-to d left-parenthesis bold x comma bold x 0 right-parenthesis plus upper O left-parenthesis StartAbsoluteValue h log h EndAbsoluteValue right-parenthesis comma

where d left-parenthesis bold x comma bold x 0 right-parenthesis is the distance function to bold x 0 .

Proof.

We prove the theorem in two dimensions. The proof can be easily extended to any number of dimensions.

First, assume the data point is a single grid point. Without loss of generality the point is at the origin. For each grid point bold x Subscript i comma j in the first quadrant its value u Subscript i comma j Superscript h only depends on its two down-left neighbors u Subscript i minus 1 comma j Superscript h Baseline comma u Subscript i comma j minus 1 Superscript h . Each grid point on the positive x axis depends on its left neighbor and each point on the positive y axis depends on its neighbor below. The first Gauss-Seidel iteration i equals 1 colon upper I comma j equals 1 colon upper J exactly propagates information from the data point to all these grid points in the right order. This order of dependence is illustrated clearly in Figure 4.1(a) for points in the first quadrant, e.g., values at grid points on the dashed line two are determined by values at grid points on dashed line one, etc. In exactly the same way grid points in the second quadrant and on the negative x axis, grid points in the third quadrant and on the negative y axis, and points in the fourth quadrant get their correct values of u Superscript h in the second, third and fourth sweep successively. Moreover the four quadrants are separated by two grid lines, i.e., the x and y axes. The values of grid points on these two lines do not depend on any values of grid points off these two lines. So the propagation of information in one quadrant during the corresponding sweep cannot cross these two lines into other quadrants. Actually, whenever a grid point achieves its correct value of u Superscript h for the system of equations (Equation2.2), it is the minimum value that can be achieved and will not change afterward.

For a data point that is not a grid point, we initially assign the exact distance values for grid points that are the vertices of the grid cell that contains the data point. The closest vertical grid line and the closest horizontal grid line partition the whole domain into four quadrants. It can be checked that the solution u Superscript h at grid points on these two lines does not depend on values of grid points off these two lines and the grid points in each quadrant get their correct values in one of the corresponding Gauss-Seidel sweeps. Figure 4.1(b) shows an example of a particular partition. This ends the proof that for a single data point the fast sweeping algorithm converges in four Gauss-Seidel iterations with alternating sweeping ordering in two dimensions.

Now we show the error estimate for grid points in the first quadrant. The proof is exactly the same for all other quadrants. Again we first assume the data point is a grid point. The exact distance function d left-parenthesis bold x right-parenthesis satisfies the Eikonal equation StartAbsoluteValue nabla d left-parenthesis bold x right-parenthesis EndAbsoluteValue equals 1 everywhere except at the data point. Using Taylor expansion at grid point bold x Subscript i comma j , we have

StartLayout 1st Row with Label left-parenthesis 4.1 right-parenthesis EndLabel StartLayout 1st Row 1st Column Blank 2nd Column d Subscript i comma j Baseline minus d Subscript i minus 1 comma j Baseline equals h left-parenthesis d Subscript x Baseline right-parenthesis Subscript i comma j Baseline minus StartFraction h squared Over 2 EndFraction d Subscript x x Baseline left-parenthesis xi Subscript i comma j Baseline right-parenthesis comma 2nd Row 1st Column Blank 2nd Column d Subscript i comma j Baseline minus d Subscript i comma j minus 1 Baseline equals h left-parenthesis d Subscript y Baseline right-parenthesis Subscript i comma j Baseline minus StartFraction h squared Over 2 EndFraction d Subscript y y Baseline left-parenthesis eta Subscript i comma j Baseline right-parenthesis comma EndLayout EndLayout

where xi Subscript i comma j and eta Subscript i comma j are two intermediate points on the line segments connecting bold x Subscript i comma j , bold x Subscript i minus 1 comma j and connecting bold x Subscript i comma j , bold x Subscript i comma j minus 1 , respectively. At bold x equals left-parenthesis x comma y right-parenthesis ,

StartLayout 1st Row with Label left-parenthesis 4.2 right-parenthesis EndLabel StartLayout 1st Row d Subscript x Baseline left-parenthesis x comma y right-parenthesis equals StartFraction x Over StartRoot x squared plus y squared EndRoot EndFraction greater-than 0 comma 0 less-than d Subscript x x Baseline left-parenthesis x comma y right-parenthesis equals StartFraction y squared Over left-parenthesis x squared plus y squared right-parenthesis Superscript 3 slash 2 Baseline EndFraction less-than-or-equal-to StartFraction 1 Over StartRoot x squared plus y squared EndRoot EndFraction comma 2nd Row d Subscript y Baseline left-parenthesis x comma y right-parenthesis equals StartFraction y Over StartRoot x squared plus y squared EndRoot EndFraction greater-than 0 comma 0 less-than d Subscript y y Baseline left-parenthesis x comma y right-parenthesis equals StartFraction x squared Over left-parenthesis x squared plus y squared right-parenthesis Superscript 3 slash 2 Baseline EndFraction less-than-or-equal-to StartFraction 1 Over StartRoot x squared plus y squared EndRoot EndFraction period EndLayout EndLayout

So the distance function satisfies the following equation at bold x Subscript i comma j :

StartLayout 1st Row with Label left-parenthesis 4.3 right-parenthesis EndLabel left-bracket d Subscript i comma j Baseline minus left-parenthesis d Subscript i minus 1 comma j Baseline minus StartFraction h squared Over 2 EndFraction d Subscript x x Baseline left-parenthesis xi Subscript i comma j Baseline right-parenthesis right-parenthesis right-bracket squared plus left-bracket d Subscript i comma j Baseline minus left-parenthesis d Subscript i comma j minus 1 Baseline minus StartFraction h squared Over 2 EndFraction d Subscript y y Baseline left-parenthesis eta Subscript i comma j Baseline right-parenthesis right-parenthesis right-bracket squared equals h squared period EndLayout

Since the solution of the quadratic equation (Equation2.3) depends monotonically on left-parenthesis a comma b right-parenthesis , we have d left-parenthesis bold x comma bold x 0 right-parenthesis less-than-or-equal-to u Superscript h Baseline left-parenthesis bold x comma bold x 0 right-parenthesis . Using the explicit expression (Equation4.2) for the derivatives and the maximum change principle, the local truncation error satisfies

StartLayout 1st Row with Label left-parenthesis 4.4 right-parenthesis EndLabel e Subscript upper T Baseline left-parenthesis bold x Subscript i comma j Baseline right-parenthesis less-than-or-equal-to StartFraction h squared Over 2 EndFraction max left-parenthesis d Subscript x x Baseline left-parenthesis bold x Subscript i minus 1 comma j Baseline right-parenthesis comma d Subscript y y Baseline left-parenthesis bold x Subscript i comma j minus 1 Baseline right-parenthesis right-parenthesis less-than-or-equal-to StartFraction h squared Over 2 min left-parenthesis d left-parenthesis bold x Subscript i minus 1 comma j Baseline right-parenthesis comma d left-parenthesis bold x Subscript i comma j minus 1 Baseline right-parenthesis right-parenthesis EndFraction period EndLayout

The global error estimate comes from the fact the accumulation of truncation errors is in the same direction of information propagation as is shown in Figure 4.1(a); i.e., grid points on line i plus j equals k depend only on grid points on line i plus j less-than-or-equal-to k minus 1 . Define e Subscript k Baseline equals max Underscript i plus j equals k Endscripts left-parenthesis u Subscript i comma j Superscript h Baseline minus d Subscript i comma j Baseline right-parenthesis . Using the simple fact that

d left-parenthesis bold x Subscript i comma j Baseline right-parenthesis equals h StartRoot i squared plus j squared EndRoot comma StartFraction left-parenthesis i plus j right-parenthesis squared Over 2 EndFraction less-than-or-equal-to i squared plus j squared less-than-or-equal-to left-parenthesis i plus j right-parenthesis squared comma

and the maximum change principle, we have

StartLayout 1st Row 1st Column u Subscript i comma j Superscript h Baseline minus d Subscript i comma j 2nd Column less-than-or-equal-to max left-parenthesis u Subscript i minus 1 comma j Superscript h Baseline minus d Subscript i minus 1 comma j Baseline plus StartFraction h squared Over 2 EndFraction d Subscript x x Baseline left-parenthesis bold x Subscript i minus 1 comma j Baseline right-parenthesis comma 2nd Row 1st Column Blank 2nd Column u Subscript i comma j minus 1 Superscript h Baseline minus d Subscript i comma j minus 1 Baseline plus StartFraction h squared Over 2 EndFraction d Subscript x x Baseline left-parenthesis bold x Subscript i comma j minus 1 Baseline right-parenthesis right-parenthesis 3rd Row 1st Column Blank 2nd Column less-than-or-equal-to e Subscript k minus 1 Baseline plus StartFraction h squared Over 2 min left-parenthesis d left-parenthesis bold x Subscript i minus 1 comma j Baseline right-parenthesis comma d left-parenthesis bold x Subscript i comma j minus 1 Baseline right-parenthesis right-parenthesis EndFraction less-than-or-equal-to e Subscript k minus 1 Baseline plus StartFraction h Over StartRoot 2 EndRoot left-parenthesis k minus 1 right-parenthesis EndFraction period EndLayout

So the maximum error that can be accumulated at bold x Subscript i comma j for k equals i plus j is

StartLayout 1st Row with Label left-parenthesis 4.5 right-parenthesis EndLabel e Subscript k Baseline less-than-or-equal-to e 1 plus StartFraction h Over StartRoot 2 EndRoot EndFraction sigma-summation Underscript k equals 1 Overscript i plus j minus 1 Endscripts StartFraction 1 Over k EndFraction less-than-or-equal-to e 1 plus StartFraction h Over StartRoot 2 EndRoot EndFraction left-parenthesis 1 plus ln left-parenthesis i plus j minus 1 right-parenthesis right-parenthesis equals upper O left-parenthesis h log left-parenthesis StartFraction 1 Over h EndFraction right-parenthesis right-parenthesis period EndLayout

Here e 1 is the maximum error for grid points on the line i plus j equals 1 , which is 0 if the data point is a grid point. The proof and estimate are exactly the same in n dimensions.

If the data point bold x 0 is not a grid point, all the error estimates are the same except that e 1 equals upper O left-parenthesis h right-parenthesis , which yields the same result.

Remark

The error estimate is sharp since there is no cancellation of local truncation errors and the accumulation of truncation errors for grid points bold x Subscript i comma i on the diagonal is exactly of upper O left-parenthesis h log left-parenthesis StartFraction 1 Over h EndFraction right-parenthesis right-parenthesis .

When there is more than one data point, the situation becomes more complicated because there are interactions among data points. Characteristics from different data points intersect and the distance function is not differentiable at equal distance points. For the exact distance function to a data set composed of discrete points, i.e., normal upper Gamma equals StartSet bold x Subscript m Baseline EndSet Subscript m equals 1 Superscript upper M , the interaction is simply the minimum rule, i.e.,

d left-parenthesis bold x comma normal upper Gamma right-parenthesis equals min left-bracket d left-parenthesis bold x comma bold x 1 right-parenthesis comma d left-parenthesis bold x comma bold x 2 right-parenthesis comma period period period comma d left-parenthesis bold x comma bold x Subscript upper M Baseline right-parenthesis right-bracket period

Let u Superscript h Baseline left-parenthesis bold x comma bold x Subscript i Baseline right-parenthesis be the numerical solution to the distance function to a single point bold x Subscript i by the fast sweeping method and define

StartLayout 1st Row with Label left-parenthesis 4.6 right-parenthesis EndLabel ModifyingBelow u With bar Superscript h Baseline left-parenthesis bold x comma normal upper Gamma right-parenthesis equals min left-bracket u Superscript h Baseline left-parenthesis bold x comma bold x 1 right-parenthesis comma u Superscript h Baseline left-parenthesis bold x comma bold x 2 right-parenthesis comma period period period comma u Superscript h Baseline left-parenthesis bold x comma bold x Subscript upper M Baseline right-parenthesis right-bracket period EndLayout

From our previous results for a single point we have

d left-parenthesis bold x comma normal upper Gamma right-parenthesis less-than-or-equal-to ModifyingBelow u With bar Superscript h Baseline left-parenthesis bold x comma normal upper Gamma right-parenthesis less-than-or-equal-to d left-parenthesis bold x comma normal upper Gamma right-parenthesis plus upper O left-parenthesis StartAbsoluteValue h log h EndAbsoluteValue right-parenthesis

after four sweeps.

Lemma 4.2

For an arbitrary set of discrete points normal upper Gamma equals StartSet bold x Subscript m Baseline EndSet Subscript m equals 1 Superscript upper M , u Superscript h Baseline left-parenthesis bold x comma normal upper Gamma right-parenthesis less-than-or-equal-to ModifyingBelow u With bar Superscript h Baseline left-parenthesis bold x comma normal upper Gamma right-parenthesis .

Proof.

For any fixed bold x there is an i , 1 less-than-or-equal-to i less-than-or-equal-to upper M , such that

ModifyingBelow u With bar Superscript h Baseline left-parenthesis bold x comma normal upper Gamma right-parenthesis equals min left-bracket u Superscript h Baseline left-parenthesis bold x comma bold x 1 right-parenthesis comma u Superscript h Baseline left-parenthesis bold x comma bold x 2 right-parenthesis comma period period period comma u Superscript h Baseline left-parenthesis bold x comma bold x Subscript upper M Baseline right-parenthesis right-bracket equals u Superscript h Baseline left-parenthesis bold x comma bold x Subscript i Baseline right-parenthesis period

After the initialization step, u Superscript h Baseline left-parenthesis bold x comma normal upper Gamma right-parenthesis less-than-or-equal-to u Superscript h Baseline left-parenthesis bold x comma bold x Subscript i Baseline right-parenthesis comma 1 less-than-or-equal-to i less-than-or-equal-to upper M . From the monotonicity in initial data for the fast sweeping algorithm, stated in Lemma 3.3, we have

u Superscript h Baseline left-parenthesis bold x comma normal upper Gamma right-parenthesis less-than-or-equal-to u Superscript h Baseline left-parenthesis bold x comma bold x Subscript i Baseline right-parenthesis equals ModifyingBelow u With bar Superscript h Baseline left-parenthesis bold x comma normal upper Gamma right-parenthesis

after any number of sweeps.

Let ModifyingAbove u With bar Superscript h Baseline left-parenthesis bold x comma normal upper Gamma right-parenthesis be the solution to the system of discretized equations, e.g., (Equation2.2) in two dimensions.

Theorem 4.3

For an arbitrary set of discrete points normal upper Gamma equals StartSet bold x Subscript m Baseline EndSet Subscript m equals 1 Superscript upper M , the numerical solution u Superscript h Baseline left-parenthesis bold x comma normal upper Gamma right-parenthesis by the fast sweeping method after 2 Superscript n sweeps, satisfies

ModifyingAbove u With bar Superscript h Baseline left-parenthesis bold x comma normal upper Gamma right-parenthesis less-than-or-equal-to u Superscript h Baseline left-parenthesis bold x comma normal upper Gamma right-parenthesis less-than-or-equal-to d left-parenthesis bold x comma normal upper Gamma right-parenthesis plus upper O left-parenthesis StartAbsoluteValue h log h EndAbsoluteValue right-parenthesis period

Proof.

The solution to the system of discretized equations, ModifyingAbove u With bar Superscript h Baseline left-parenthesis bold x comma normal upper Gamma right-parenthesis , can be viewed as the solution by the fast sweeping algorithm after the iteration converges as is shown in Theorem 3.6. Since the solution of the fast sweeping algorithm is nonincreasing with Gauss-Seidel iterations, we have u Superscript h Baseline left-parenthesis bold x comma normal upper Gamma right-parenthesis greater-than-or-equal-to ModifyingAbove u With bar Superscript h Baseline left-parenthesis bold x comma normal upper Gamma right-parenthesis after any number of sweeps. So after 2 Superscript n sweeps, the numerical solution u Superscript h Baseline left-parenthesis bold x comma normal upper Gamma right-parenthesis produced by the fast sweeping algorithm satisfies

ModifyingAbove u With bar Superscript h Baseline left-parenthesis bold x comma normal upper Gamma right-parenthesis less-than-or-equal-to u Superscript h Baseline left-parenthesis bold x comma normal upper Gamma right-parenthesis less-than-or-equal-to ModifyingBelow u With bar Superscript h Baseline left-parenthesis bold x comma normal upper Gamma right-parenthesis equals d left-parenthesis bold x comma normal upper Gamma right-parenthesis plus upper O left-parenthesis h log StartFraction 1 Over h EndFraction right-parenthesis period

Since the upwind difference is of first order accuracy, StartAbsoluteValue ModifyingAbove u With bar Superscript h Baseline left-parenthesis bold x comma normal upper Gamma right-parenthesis minus d left-parenthesis bold x comma normal upper Gamma right-parenthesis EndAbsoluteValue is at most upper O left-parenthesis h right-parenthesis . The general results for Hamilton Jacobi equations, e.g., Reference5Reference14Reference3Reference7, show that the numerical solution from a consistent and monotone scheme converges to the viscosity solution with the order of h Superscript one-half . The upper bound in the above theorem is sharp as is shown in Theorem 4.1. If the error estimate, StartAbsoluteValue ModifyingAbove u With bar Superscript h Baseline left-parenthesis bold x comma normal upper Gamma right-parenthesis minus d left-parenthesis bold x comma normal upper Gamma right-parenthesis EndAbsoluteValue , is also upper O left-parenthesis StartAbsoluteValue h log h EndAbsoluteValue right-parenthesis or worse, the theorem says that for the distance function, the iterative solution after 2 Superscript n sweeps is as accurate as ModifyingAbove u With bar Superscript h Baseline left-parenthesis bold x comma normal upper Gamma right-parenthesis . Any other method that solves the same discretized system of equations has the same accuracy too.

However, we do not have d left-parenthesis bold x comma normal upper Gamma right-parenthesis less-than-or-equal-to u Superscript h Baseline left-parenthesis bold x comma normal upper Gamma right-parenthesis for a general data set normal upper Gamma due to the interactions among data points. Figure 4.2 shows an example of two data points. At those circled grid points the characteristics from both data points meet. The distance function follows either one of the characteristics. For example at grid point (1, 1) the distance function satisfies

left-bracket left-parenthesis d Subscript 1 comma 1 Baseline minus d Subscript 1 comma 2 Baseline right-parenthesis Superscript plus Baseline right-bracket squared equals h squared or left-bracket left-parenthesis d Subscript 1 comma 1 Baseline minus d Subscript 2 comma 1 Baseline right-parenthesis Superscript plus Baseline right-bracket squared equals h squared

and d Subscript 1 comma 1 Baseline equals 2 h . However in our upwind scheme, the numerical solution u Superscript h uses both characteristics and satisfies

left-bracket left-parenthesis u Subscript 1 comma 1 Superscript h Baseline minus u Subscript 1 comma 2 Superscript h Baseline right-parenthesis Superscript plus Baseline right-bracket squared plus left-bracket left-parenthesis u Subscript 1 comma 1 Superscript h Baseline minus u Subscript 2 comma 1 Superscript h Baseline right-parenthesis Superscript plus Baseline right-bracket squared equals h squared comma

which gives u Subscript 1 comma 1 Superscript h Baseline equals left-parenthesis 1 plus StartFraction 1 Over StartRoot 2 EndRoot EndFraction right-parenthesis h less-than d Subscript 1 comma 1 . We can view this as the information propagation speed is numerically doubled. However, since the upwind scheme uses at most two characteristics from u Subscript x min in the x -direction and from u Subscript y min in the y -direction in two dimensions, we show that this is actually the worst truncation error that can occur at a grid point due to the interactions of data points, i.e. when characteristics intersect orthogonally and align with both axes. For instance in two dimensions, without loss of generality suppose u Subscript i comma j Superscript h Baseline greater-than-or-equal-to u Subscript i comma j minus 1 Superscript h Baseline greater-than-or-equal-to u Subscript i minus 1 comma j Superscript h and

left-parenthesis u Subscript i comma j Superscript h Baseline minus u Subscript i minus 1 comma j Superscript h Baseline right-parenthesis squared plus left-parenthesis u Subscript i comma j Superscript h Baseline minus u Subscript i comma j minus 1 Superscript h Baseline right-parenthesis squared equals h squared period

Then left-parenthesis u Subscript i comma j Superscript h Baseline minus u Subscript i minus 1 comma j Superscript h Baseline right-parenthesis squared greater-than-or-equal-to h squared slash 2 and the equality holds when u Subscript i comma j minus 1 Superscript h Baseline equals u Subscript i minus 1 comma j Superscript h . On the other hand the distance function satisfies left-parenthesis d Subscript i comma j Baseline minus d Subscript i minus 1 comma j Baseline right-parenthesis squared less-than-or-equal-to h squared and the equality holds when the x axis is a characteristic. In n dimensions, the characteristics can be used at most n times when they intersect at one grid point. So the worst local truncation error due to the interactions of characteristics is StartRoot 1 minus StartFraction 1 Over n EndFraction EndRoot h . For the modified version of the fast sweeping algorithm (Equation2.7), the truncation errors at equal distance points can be twice as much.

To get a clearer picture of the convergence of the iteration and error estimate for a general data set, we have to study the interactions among data points more carefully. We can partition all grid points into the Voronoi cell of each data point. The Voronoi diagram is according to the the numerical solution u Superscript h Baseline left-parenthesis bold x comma bold x 1 right-parenthesis comma u Superscript h Baseline left-parenthesis bold x comma bold x 2 right-parenthesis comma period period period comma u Superscript h Baseline left-parenthesis bold x comma bold x Subscript upper M Baseline right-parenthesis ; i.e., a grid point bold x is in the Voronoi cell of bold x Subscript m if

u Superscript h Baseline left-parenthesis bold x comma bold x Subscript m Baseline right-parenthesis equals min left-bracket u Superscript h Baseline left-parenthesis bold x comma bold x 1 right-parenthesis comma u Superscript h Baseline left-parenthesis bold x comma bold x 2 right-parenthesis comma period period period comma u Superscript h Baseline left-parenthesis bold x comma bold x Subscript upper M Baseline right-parenthesis right-bracket period

If a grid point and all its neighboring grid points belong to the same Voronoi cell, we call it an interior point. Otherwise we call it a boundary point. The interaction of different data points occurs only at boundary points. Figure 4.3(a) shows a typical Voronoi cell for a data point bold x Subscript m . For cell boundary points (those circled points), u Superscript h Baseline left-parenthesis bold x comma normal upper Gamma right-parenthesis may pick up information from more than one data point.

To get the lower bound for the numerical solution after 2 Superscript n sweeps, we use ModifyingBelow u With bar Superscript h Baseline left-parenthesis bold x comma normal upper Gamma right-parenthesis equals min left-bracket u Superscript h Baseline left-parenthesis bold x comma bold x 1 right-parenthesis comma u Superscript h Baseline left-parenthesis bold x comma bold x 2 right-parenthesis comma period period period comma u Superscript h Baseline left-parenthesis bold x comma bold x Subscript upper M Baseline right-parenthesis right-bracket as the initial data and start the fast sweeping iteration. Due to the monotonicity in initial data we get a solution that provides a lower bound for the numerical solution for which we use the standard initialization step. ModifyingBelow u With bar Superscript h Baseline left-parenthesis bold x comma normal upper Gamma right-parenthesis already satisfies the discretized equations at interior points of each Voronoi cell. After we start the fast sweeping algorithm, the decrease of the values at interior points of each Voronoi cell is caused by the interactions at Voronoi cell boundaries. Moreover, if we start with ModifyingBelow u With bar Superscript h Baseline left-parenthesis bold x Subscript i comma j Baseline comma normal upper Gamma right-parenthesis , it is easy to show StartAbsoluteValue ModifyingBelow u With bar Superscript h Baseline left-parenthesis bold x Subscript i comma j Baseline comma normal upper Gamma right-parenthesis minus ModifyingBelow u With bar Superscript h Baseline left-parenthesis bold x Subscript i plus-or-minus 1 comma j plus-or-minus 1 Baseline comma normal upper Gamma right-parenthesis EndAbsoluteValue less-than-or-equal-to h from the system of discretized equations at each grid point and the definition of Voronoi cells. Hence from the maximum change principle Proposition 3.2, we can imagine that the maximum decrease of values at all grid points due to the interactions at the Voronoi cell boundary is of order h . But unlike the case for the real distance function where information propagates only along characteristics and all characteristics flow into the Voronoi cell boundary, in the finite difference scheme a grid point may have a larger domain of dependence as is illustrated in Figure 4.3(b). So interactions at Voronoi cell boundaries may propagate into the cell. This may also cause more than 2 Superscript n sweeps for convergence in n dimensions. Example 1 in Section 6 shows that even for two data points in two dimensions, more than four sweeps are needed for the iteration to converge.

Now we consider computing the distance function to an arbitrary set. For example, instead of discrete points, normal upper Gamma is a smooth curve or surface.

Theorem 4.4

If the distance function in the neighborhood of an arbitrary data set normal upper Gamma in upper R Superscript n is given initially, let u Superscript h Baseline left-parenthesis bold x comma normal upper Gamma right-parenthesis be the numerical solution by the fast sweeping method after 2 Superscript n sweeps. We have

ModifyingAbove u With bar Superscript h Baseline left-parenthesis bold x comma normal upper Gamma right-parenthesis less-than-or-equal-to u Superscript h Baseline left-parenthesis bold x comma normal upper Gamma right-parenthesis less-than-or-equal-to d left-parenthesis bold x comma normal upper Gamma right-parenthesis plus upper O left-parenthesis h log StartFraction 1 Over h EndFraction right-parenthesis comma

where ModifyingAbove u With bar Superscript h Baseline left-parenthesis bold x comma normal upper Gamma right-parenthesis is the solution to the discretized system Equation2.2.

Proof.

Let normal upper Gamma overbar be the set of grid points that encloses the set normal upper Gamma ; i.e., normal upper Gamma overbar contains vertices of all those grid cells that intersect with normal upper Gamma . We have

StartLayout 1st Row with Label left-parenthesis 4.7 right-parenthesis EndLabel StartAbsoluteValue d left-parenthesis bold x comma normal upper Gamma overbar right-parenthesis minus d left-parenthesis bold x comma normal upper Gamma right-parenthesis EndAbsoluteValue equals upper O left-parenthesis h right-parenthesis comma EndLayout

since for any bold y element-of normal upper Gamma comma there-exists bold y Subscript i comma j Baseline element-of normal upper Gamma overbar such that StartAbsoluteValue bold y Subscript i comma j Baseline minus bold y EndAbsoluteValue equals upper O left-parenthesis h right-parenthesis and vice versa.

By the monotonicity in initial data, u Superscript h Baseline left-parenthesis bold x comma normal upper Gamma right-parenthesis greater-than-or-equal-to u Superscript h Baseline left-parenthesis bold x comma normal upper Gamma overbar right-parenthesis after any number of sweeps, since initially u Superscript h Baseline left-parenthesis bold x comma normal upper Gamma right-parenthesis starts with distance to normal upper Gamma for bold x element-of normal upper Gamma overbar while u Superscript h Baseline left-parenthesis bold x comma normal upper Gamma overbar right-parenthesis equals 0 for bold x element-of normal upper Gamma overbar . However, the initial difference between u Superscript h Baseline left-parenthesis bold x comma normal upper Gamma right-parenthesis and u Superscript h Baseline left-parenthesis bold x comma normal upper Gamma overbar right-parenthesis is upper O left-parenthesis h right-parenthesis . By the contraction property from Corollary 3.5, we have

StartLayout 1st Row with Label left-parenthesis 4.8 right-parenthesis EndLabel u Superscript h Baseline left-parenthesis bold x comma normal upper Gamma right-parenthesis minus u Superscript h Baseline left-parenthesis bold x comma normal upper Gamma overbar right-parenthesis equals upper O left-parenthesis h right-parenthesis comma for-all bold x comma EndLayout

after any number of sweeps. We apply Theorem 4.3 to u Superscript h Baseline left-parenthesis bold x comma normal upper Gamma overbar right-parenthesis and combine it with (Equation4.7) and (Equation4.8) to finish the proof.

Actually for an arbitrary data set, which can be discrete points and/or continuous manifolds, we only need approximate distance values at grid points near the data set within first order accuracy, since the upwind finite difference scheme is at most of first order.

5. General Eikonal equations

For the general Eikonal equation (Equation2.1), the characteristics are curves starting from the boundary. The key issue is the maximum number of sweeps needed to cover information propagation along a single characteristic curve. This number, which is analogous to the condition number for elliptic equations, determines the number of iterations needed for the fast sweeping algorithm. For a single characteristic curve starting at a point bold x 0 element-of normal upper Gamma , we divide it into the least number of pieces such that the tangent directions in each piece belong to the same quadrant. Information propagation along each piece can be covered by one of the sweepings ordering successively. So the number of sweeps needed to cover the whole characteristic curve is proportional to the number of pieces or turns. Figure 5.1 shows an example in 2D. The characteristic curve starting from bold x 0 element-of normal upper Gamma can be divided into five pieces. The tangent directions in each piece belong to the same quadrant. If we order one round of four alternating sweeps as in Section 2, the first and the second pieces are covered by the first and fourth sweeps in the first round, respectively, the third piece is covered by the third sweep in the second round, the fourth piece is covered by the second sweep in the third round and the fifth piece is covered by the first sweep in the fourth round.

One quantity that can characterize how sharp the tangent of a curve can turn is curvature. The following lemma shows a bound on the curvature of any characteristic curve.

Lemma 5.1

The maximum curvature for any characteristic curve of equation Equation2.1 is bounded by max Underscript bold x element-of normal upper Omega Endscripts StartAbsoluteValue StartFraction nabla f left-parenthesis bold x right-parenthesis Over f left-parenthesis bold x right-parenthesis EndFraction EndAbsoluteValue .

Proof.

Denote upper H left-parenthesis bold q comma bold x right-parenthesis equals StartAbsoluteValue bold q EndAbsoluteValue minus f left-parenthesis bold x right-parenthesis , where bold q equals nabla u . The characteristic equation is

StartLayout Enlarged left-brace 1st Row ModifyingAbove bold x With dot equals nabla Subscript bold q Baseline upper H equals StartFraction nabla u Over f left-parenthesis bold x right-parenthesis EndFraction comma 2nd Row ModifyingAbove bold q With dot equals minus nabla Subscript bold x Baseline upper H equals nabla f left-parenthesis bold x right-parenthesis comma 3rd Row ModifyingAbove u With dot equals nabla u dot ModifyingAbove bold x With dot equals f left-parenthesis bold x right-parenthesis period EndLayout

The information propagates along the characteristics from smaller u to larger u . Since StartAbsoluteValue ModifyingAbove bold x With dot EndAbsoluteValue equals 1 , the curvature along a characteristic is

StartLayout 1st Row 1st Column ModifyingAbove bold x With two-dots 2nd Column equals StartFraction nabla ModifyingAbove u With dot Over f left-parenthesis bold x right-parenthesis EndFraction minus StartFraction nabla u Over f squared left-parenthesis bold x right-parenthesis EndFraction nabla f dot ModifyingAbove bold x With dot 2nd Row 1st Column Blank 2nd Column equals StartFraction nabla f left-parenthesis bold x right-parenthesis Over f left-parenthesis bold x right-parenthesis EndFraction minus left-parenthesis StartFraction nabla f left-parenthesis bold x right-parenthesis Over f left-parenthesis bold x right-parenthesis EndFraction dot StartFraction nabla u Over StartAbsoluteValue nabla u EndAbsoluteValue EndFraction right-parenthesis StartFraction nabla u Over StartAbsoluteValue nabla u EndAbsoluteValue EndFraction 3rd Row 1st Column Blank 2nd Column equals StartFraction nabla f left-parenthesis bold x right-parenthesis Over f left-parenthesis bold x right-parenthesis EndFraction minus upper P Subscript bold n Baseline StartFraction nabla f left-parenthesis bold x right-parenthesis Over f left-parenthesis bold x right-parenthesis EndFraction comma EndLayout

where upper P Subscript bold n is the projection on the normal direction bold n equals StartFraction nabla u Over StartAbsoluteValue nabla u EndAbsoluteValue EndFraction . So StartAbsoluteValue ModifyingAbove bold x With two-dots EndAbsoluteValue less-than-or-equal-to StartAbsoluteValue StartFraction nabla f left-parenthesis bold x right-parenthesis Over f left-parenthesis bold x right-parenthesis EndFraction EndAbsoluteValue .

So for general Eikonal equations the number of iterations for the fast sweeping method depends on the right-hand side f left-parenthesis bold x right-parenthesis and the size and dimension of the computational domain only. The computed numerical solution has the same accuracy as the solution by any other method that uses the same discretization.

If the boundary normal upper Gamma and f left-parenthesis bold x right-parenthesis are smooth and f left-parenthesis bold x right-parenthesis greater-than 0 , then normal upper Gamma is a noncharacteristic boundary and there is a neighborhood of normal upper Gamma in which characteristics do not cross each other and the solution u left-parenthesis bold x right-parenthesis is smooth (see Reference8). Let this neighborhood be normal upper Omega Subscript normal upper Gamma . We have

Theorem 5.2

The numerical solution u Subscript i comma j Superscript h to the discretized system Equation2.2 is of first order in normal upper Omega Subscript normal upper Gamma ; i.e., StartAbsoluteValue u Subscript i comma j Superscript h Baseline minus u left-parenthesis bold x Subscript i comma j Baseline right-parenthesis EndAbsoluteValue equals upper O left-parenthesis h right-parenthesis comma bold x Subscript i comma j Baseline element-of normal upper Omega Subscript normal upper Gamma Baseline .

Proof.

Without loss of generality, suppose the numerical solution to (Equation2.2) satisfies the equation

left-parenthesis u Subscript i comma j Superscript h Baseline minus u Subscript i minus 1 comma j Superscript h Baseline right-parenthesis squared plus left-parenthesis u Subscript i comma j Superscript h Baseline minus u Subscript i comma j minus 1 Superscript h Baseline right-parenthesis squared equals f Subscript i comma j Superscript 2 Baseline h squared

at a grid point bold x Subscript i comma j , while the true solution u left-parenthesis bold x right-parenthesis satisfies

left-bracket u Subscript i comma j Baseline minus left-parenthesis u Subscript i minus 1 comma j Baseline minus StartFraction h squared Over 2 EndFraction u Subscript x x Baseline left-parenthesis xi Subscript i comma j Baseline right-parenthesis right-parenthesis right-bracket squared plus left-bracket u Subscript i comma j Baseline minus left-parenthesis u Subscript i comma j minus 1 Baseline minus StartFraction h squared Over 2 EndFraction u Subscript y y Baseline left-parenthesis eta Subscript i comma j Baseline right-parenthesis right-parenthesis right-bracket squared equals f Subscript i comma j Superscript 2 Baseline h squared

at bold x Subscript i comma j , where xi Subscript i comma j and eta Subscript i comma j are two intermediate points on the line segments connecting bold x Subscript i comma j , bold x Subscript i minus 1 comma j and connecting bold x Subscript i comma j , bold x Subscript i comma j minus 1 , respectively. Since u Subscript x x Baseline comma u Subscript y y Baseline are bounded, from the maximum change property, Lemma 3.1, we can deduce that the local truncation error is upper O left-parenthesis h squared right-parenthesis . The propagation and accumulation of truncation errors following the causality along characteristics in a finite domain is at most upper O left-parenthesis h right-parenthesis .

This error estimate breaks down when the solution has singularities. Since the characteristics do intersect for general Eikonal equations, we cannot use this argument after characteristics intersect. We quote the general error estimate results for monotone schemes for Hamilton-Jacobi equations Reference5Reference14Reference3Reference7. The error is of order upper O left-parenthesis h Superscript 1 slash 2 Baseline right-parenthesis . In general, there are two scenarios for the solution to have singularities. In the first scenario normal upper Gamma is smooth but characteristics intersect and shocks are formed. The solution is continuous but not differentiable at shocks. The numerical solution can be only first order accurate at shocks. However, characteristics and information flow into shocks for the true solution. Numerically we also observe that errors made at shocks do not propagate away from shocks and hence high order schemes can achieve high order accuracy in smooth regions. In the second scenario normal upper Gamma has singularities such as corners and kinks. Hence the solution also has singularities at normal upper Gamma . The distance function to a single point is such an example. Again only first order accuracy can be achieved near the boundary for any numerical scheme. Since characteristics and information flow out of the boundary, errors made near the boundary will propagate out to the computational domain. So the global error will be at most of first order no matter what scheme is used. The only solution is to use a finer grid near singularities at the boundary.

6. Numerical results

In this section we will use numerical examples to test the fast sweeping algorithm and to verify the analysis in previous sections. We can compute the distance function only in a narrow neighborhood of the data set as is described in Section 2 if needed, which saves an order of magnitude of computational cost.

For computing the distance function to a set of discrete points, we use the following procedure for the initialization step. First initialize the distance value of all grid points to be a large value, which should be larger than the maximum possible values for our later computed distance value u Superscript h in the domain. Then go through each data point and update the distance values of its neighboring grid points. For example, for each data point we find the grid cell that contains it and then compute the exact distance value of vertices of the grid cell to the data point. We replace the current values of these vertices whenever distance to this data point provides a smaller value. Of course we can include more neighboring grid points for which the exact distance values are computed. In our calculations, distance values are computed at grid points that are within two grid cells of the data set in the initialization step except for the first example. After going through all data points, we have computed exact distance values at those grid points in a neighborhood of the data set. All other grid points remain to have a large value. This procedure is of complexity upper O left-parenthesis upper M right-parenthesis for upper M data points. In general we can find the global distance function to any data set as long as the distance values on grid points neighboring the set are provided or computed initially.

Example 1

This example shows the interaction of two data points on a simple grid in two dimensions. Five iterations are needed for the fast sweeping algorithm to converge. However, changes after four iterations are upper O left-parenthesis h right-parenthesis no matter what size the grid is as we have tested. In Figure 6.1 we show the results after each iteration on a 7 times 7 grid. The two data points are grid points at left-parenthesis 2 comma 6 right-parenthesis comma left-parenthesis 5 comma 2 right-parenthesis . For this example, we scale the grid size to be 1. Initially, as is shown in Table 6.1(a), we assign a large enough number (100 is enough for this grid) to grid points that are not data points and assign zero to those two data points. Table 6.1(b) shows the numerical solution after the first sweep, i equals 1 colon 7 comma j equals 1 colon 7 . Table 6.1(c)–(f) shows the numerical solution after second, third, fourth and fifth sweep. Table 6.2 shows the maximum change between each sweep. The changes in the first four sweeps are significant since every time there are grid points whose values change from their initial (large) assignments to the correct values. The change between fourth sweep and fifth sweep, which is caused by the interaction of two points when their characteristics intersect at the Voronoi cell boundary, is much less than upper O left-parenthesis h right-parenthesis ( h equals 1 in this case). For tests with different grid sizes and different locations of two data points, it may take more than five iterations to converge but changes after four iterations are always small. Table 6.1(g) shows the exact distance function. Those underlined numerical values in Table 6.1(e) show where the numerical solution is smaller than the exact distance function due to the interaction of these two data points. Table 6.1(h) shows the Voronoi diagram according to the numerical solutions of the distance function to each single data point, as is explained in Section 4. The integer at each grid point shows to which data point it is closer. We see the interaction occurs exactly at the Voronoi boundary. All these numerical results agree with our analysis.

Example 2

In this example we compute the distance function to discrete points in both two and three dimensions to test the convergence and accuracy of the fast sweeping method. In the first case we have a single data point in both two and three dimensions. The domain is a unit square/cube. The data point is located at left-parenthesis StartFraction 1 Over StartRoot 2 EndRoot EndFraction comma StartFraction 1 Over StartRoot 3 EndRoot EndFraction right-parenthesis in two dimensions and left-parenthesis StartFraction 1 Over StartRoot 2 EndRoot EndFraction comma StartFraction 1 Over StartRoot 3 EndRoot EndFraction comma 0.1 pi right-parenthesis in three dimensions. Table 6.3 shows errors measured in different norms with different grid sizes. For one single point the fast sweeping method converges exactly in four sweeps in two dimensions and in eight sweeps in three dimensions.

In the second case, we generate a set of 100 random points in a unit square in two dimensions and in a unit cube in three dimensions. In two dimensions it takes up to eight sweeps to converge. In three dimensions it takes up to 20 sweeps to converge. We show in Table 6.4 the errors after four sweeps in two dimensions and eight sweeps in three dimensions.

Example 3

In this example we compute the distance function to two continuous sets, four linked circles in two dimensions and four spheres in three dimensions. Again it takes more than four sweeps in two dimensions and eight sweeps in three dimensions to converge in both cases. We show errors after four sweeps in two dimensions and eight sweeps in three dimensions in Table 6.5 In Figure 6.1 we show the contour plot of the numerical solution in two dimensions and a particular contour in three dimensions.

Example 4

In this example we present two cases for general Eikonal equations. In the first case we show convergence and order of accuracy of the fast sweeping algorithm. The exact solution is

u left-parenthesis x comma y right-parenthesis equals StartAbsoluteValue e Superscript left-parenthesis StartRoot x squared plus y squared EndRoot minus r 0 right-parenthesis left-parenthesis a x squared plus 2 b x y plus c y squared plus d right-parenthesis Baseline minus 1 EndAbsoluteValue comma

where r 0 equals 0.2 comma a equals 0.1 comma b equals 0.5 comma c equals 1 comma d equals 0.4 . Also, StartAbsoluteValue nabla u EndAbsoluteValue is computed explicitly and is used as the given f left-parenthesis x comma y right-parenthesis . The boundary condition is u left-parenthesis x comma y right-parenthesis equals 0 at the circle StartRoot x squared plus y squared EndRoot equals r 0 . We use the exact values of u left-parenthesis x comma y right-parenthesis at grid points near the circle initially. We set the convergence criterion to be double-vertical-bar u Superscript left-parenthesis k right-parenthesis Baseline minus u Superscript left-parenthesis k minus 1 right-parenthesis Baseline double-vertical-bar Subscript normal infinity Baseline less-than epsilon equals 10 Superscript negative 6 . The number of iterations and errors in different norms are shown in Table 6.6. We see that the number of iterations is independent of grid size. The contour plot of the solution is shown in Figure 6.2.

Note that the errors in all norms show first order convergence. This is due to the fact that the boundary, i.e., the initial wave front, is smooth. The only singularity is at the center where StartAbsoluteValue nabla u EndAbsoluteValue equals 0 . However, characteristics flow into the singularity; i.e., the error at the singularity is determined by the error around it. This is contrary to the case where the boundary has singularities; e.g., the boundary is a single point or there are corners at the boundary. In this case characteristics flow out of the singularities of the boundary. Hence numerical errors made at singularities will propagate out and can accumulate, which gives a global error of order h log left-parenthesis StartFraction 1 Over h EndFraction right-parenthesis just as in the case for the distance function to a single point.

In the second case we compute the first arrival time for a point source at left-parenthesis x 0 comma y 0 right-parenthesis equals left-parenthesis 0 comma 0 right-parenthesis . The Eikonal equation is

StartAbsoluteValue nabla u EndAbsoluteValue equals f left-parenthesis x comma y right-parenthesis equals 1 plus e Superscript minus 40 left-bracket left-parenthesis x plus 0.2 right-parenthesis squared plus left-parenthesis y plus 0.2 right-parenthesis squared right-bracket Baseline minus e Superscript minus 40 left-bracket left-parenthesis x minus 0.2 right-parenthesis squared plus left-parenthesis y minus 0.2 right-parenthesis squared right-bracket Baseline comma u left-parenthesis 0 comma 0 right-parenthesis equals 0 period

So the corresponding velocity field is 1 slash f left-parenthesis x comma y right-parenthesis and the ratio between the maximum velocity and the minimum velocity is 6.667 times 10 Superscript 5 . The function f left-parenthesis x comma y right-parenthesis is shown in Figure 6.3(a). In Figure 6.3(b) the solid line is the contour plot of the arrival time and the dashed line is the contour plot of f left-parenthesis x comma y right-parenthesis . The largest difference between two consecutive iterations is shown in Table 6.7. We set the initial large value to be 10 Superscript 6 . The convergence criterion is the maximum difference between two consecutive iterations less than 10 Superscript negative 6 . Although six to seven iterations are needed for this convergence criterion for different grids, we can see clearly from this table that only the first six iterations are essential. After six iterations, the difference drops dramatically and is much smaller than the grid size. This means that it takes six iterations to cover the information propagation along all characteristic curves in the computation domain. This pattern is independent of grid size. The computing time scales linearly with the grid size. For this example, it only takes 0.2 seconds for the computation on a 500 times 500 grid using a 2.4 GHz PC.

Example 5

We present a potential application for function reconstruction. In theory, for a upper C Superscript 1 function u left-parenthesis bold x right-parenthesis comma bold x element-of upper R Superscript n , if all local minima (or maxima) of u together with its gradient StartAbsoluteValue nabla u EndAbsoluteValue are known, we can reconstruct the function u left-parenthesis bold x right-parenthesis by solving the Eikonal equation with the prescribed local minima (or maxima) as the boundary condition. Here we apply this idea to the reconstruction of a discrete function. Suppose u Subscript i comma j is a two dimensional discrete function defined at grid points left-parenthesis x Subscript i Baseline comma y Subscript j Baseline right-parenthesis . We first search for all local minima. u Subscript i comma j is a local minimum if u Subscript i comma j Baseline less-than-or-equal-to min left-brace u Subscript i plus 1 comma j Baseline comma u Subscript i minus 1 comma j Baseline comma u Subscript i comma j minus 1 Baseline comma u Subscript i comma j plus 1 Baseline right-brace . We record the location, left-parenthesis i comma j right-parenthesis , and the value u Subscript i comma j . Next we extract the gradient information at points that are not local minima. In order to construct the exact inverse process for our fast sweeping algorithm, we use upwind difference to compute the gradient at a grid point left-parenthesis x Subscript i Baseline comma y Subscript j Baseline right-parenthesis as follows. Let

StartLayout 1st Row u Subscript x min Baseline equals min left-parenthesis u Subscript i minus 1 comma j Baseline comma u Subscript i plus 1 comma j Baseline right-parenthesis comma 2nd Row u Subscript y min Baseline equals min left-parenthesis u Subscript i comma j minus 1 Baseline comma u Subscript i comma j plus 1 Baseline right-parenthesis period EndLayout

Define

StartLayout 1st Row u Subscript x Baseline equals StartLayout Enlarged left-brace 1st Row 1st Column StartFraction u Subscript i comma j Baseline minus u Subscript x min Baseline Over h EndFraction 2nd Column if u Subscript i comma j Baseline greater-than u Subscript x min Baseline comma 2nd Row 1st Column 0 2nd Column if u Subscript i comma j Baseline less-than-or-equal-to u Subscript x min Baseline comma EndLayout 2nd Row u Subscript y Baseline equals StartLayout Enlarged left-brace 1st Row 1st Column StartFraction u Subscript i comma j Baseline minus u Subscript y min Baseline Over h EndFraction 2nd Column if u Subscript i comma j Baseline greater-than u Subscript y min Baseline comma 2nd Row 1st Column 0 2nd Column if u Subscript i comma j Baseline less-than-or-equal-to u Subscript y min Baseline comma EndLayout EndLayout

and

f Subscript i comma j Baseline equals StartAbsoluteValue nabla u EndAbsoluteValue Subscript i comma j Baseline equals StartRoot u Subscript x Superscript 2 Baseline plus u Subscript y Superscript 2 Baseline EndRoot period

By enforcing the values at the minima and solving the system (Equation2.2), we can recover u Subscript i comma j to machine precision. In regions where u Subscript i comma j is constant, f Subscript i comma j Baseline equals StartAbsoluteValue nabla u EndAbsoluteValue Subscript i comma j Baseline equals 0 . Here we show the reconstruction of the peak function and the hat function in Matlab. For the peak function, there are six local minimal grid points and it takes 10 iterations to converge. The reconstruction is exact to machine precision and is shown in Figure 6.4(a). For the hat function, there are many local minima due to the valleys and the number of minima scales with the number of grid points. It takes seven iterations to converge. The reconstruction is also exact to machine precision and is shown in Figure 6.4(b).

Example 6

As the last example, we present an application of the distance function in computer visualization. In Reference20, efficient algorithms are developed to analyze and visualize large sets of unorganized points using the distance function and distance contours. In particular, an appropriate distance contour can be extracted very quickly for the visualization of the data set, which avoids sorting out complicated ordering and connections among all data points. Figure 6.5 shows the visualization of a Buddha statue using a distance contour on a 156 times 371 times 156 grid from points obtained by a laser scanner. The data set has 543,652 points and is obtained from The Stanford 3D Scanning Repository. The whole process takes a few seconds due to the fast computation of the distance function.

Acknowledgment

The author would like to thank Dr. Paul Dupuis for some interesting discussions that started the work.

Figures

Figure 2.1.

The fast sweeping algorithm in one dimension.

Figure 2.1(a).

the computed distance function after the first left to right Gauss-Seidel sweeping

Figure 2.1(b).

the computed distance function after the second right to left Gauss-Seidel sweeping

Figure 2.2.

The fast sweeping algorithm in two dimensions.

Figure 2.2(a).

the fast sweeping algorithm for a single data point

Figure 2.2(b).

the fast sweeping algorithm for a circle

Figure 4.1.
Figure 4.1(a).

order of dependence in the first quadrant

Figure 4.1(b).

quadrant partition for an arbitrary data point

Figure 4.2.

The numerical solution for two data points bold x 1 comma bold x 2

Figure 4.3.
Figure 4.3(a).

error at the cell boundary and its propagation

Figure 4.3(b).

domain of dependence of a grid point

Figure 5.1.

Division of a characteristic curve for a general Eikonal equation.

Table 6.1.
Table 6.1(a).

initial setup

Table 6.1(b).

after 1st sweep

Table 6.1(c).

after 2nd sweep

Table 6.1(d).

after 3rd sweep

Table 6.1(e).

after 4th sweep

Table 6.1(f).

after 5th sweep

Table 6.1(g).

exact distance

Table 6.1(h).

Voronoi diagram

Table 6.2.

maximum change in each sweep

Table 6.3.

distance function to a single data point

Table 6.3(a).

two dimensions

Table 6.3(b).

three dimensions

Table 6.4.

distance function to a set of 100 random data points

Table 6.4(a).

two dimensions

Table 6.4(b).

three dimensions

Table 6.5.

distance function to a continuous set

Table 6.5(a).

two dimensions

Table 6.5(b).

three dimensions

Figure 6.1.
Figure 6.1(a).

contour plot in 2D

Figure 6.1(b).

the contour u Superscript h Baseline equals 0.03 in 3D

Figure 6.2.

Contour plot of the solution; dotted line is where the boundary condition is prescribed.

Figure 6.3.
Figure 6.3(a).

graph of f left-parenthesis x comma y right-parenthesis

Figure 6.3(b).

contour plot of the solution (solid line) and f left-parenthesis x comma y right-parenthesis (dashed line)

Table 6.6.

Table 6.7.

maximum difference between two consecutive iterations

Figure 6.4.
Figure 6.4(a).

reconstruction of the peak function

Figure 6.4(b).

reconstruction of the hat function

Figure 6.5.

Visualization of 3D scanned points using a distance contour.

Figure 6.5(a).

front

Figure 6.5(b).

diagonal

Figure 6.5(c).

back

Mathematical Fragments

Equation (2.1)
StartLayout 1st Row with Label left-parenthesis 2.1 right-parenthesis EndLabel StartLayout 1st Row 1st Column StartAbsoluteValue nabla u left-parenthesis bold x right-parenthesis EndAbsoluteValue equals f left-parenthesis bold x right-parenthesis comma 2nd Column bold x element-of upper R Superscript n Baseline comma 2nd Row 1st Column u left-parenthesis bold x right-parenthesis equals 0 comma 2nd Column bold x element-of normal upper Gamma subset-of upper R Superscript n Baseline comma EndLayout EndLayout
Equation (2.2)
StartLayout 1st Row with Label left-parenthesis 2.2 right-parenthesis EndLabel StartLayout 1st Row left-bracket left-parenthesis u Subscript i comma j Superscript h Baseline minus u Subscript x min Superscript h Baseline right-parenthesis Superscript plus Baseline right-bracket squared plus left-bracket left-parenthesis u Subscript i comma j Superscript h Baseline minus u Subscript y min Superscript h Baseline right-parenthesis Superscript plus Baseline right-bracket squared equals f Subscript i comma j Superscript 2 Baseline h squared 2nd Row Blank 3rd Row i equals 2 comma period period period comma upper I minus 1 comma j equals 2 comma period period period comma upper J minus 1 comma EndLayout EndLayout
Equation (2.3)
StartLayout 1st Row with Label left-parenthesis 2.3 right-parenthesis EndLabel left-bracket left-parenthesis x minus a right-parenthesis Superscript plus Baseline right-bracket squared plus left-bracket left-parenthesis x minus b right-parenthesis Superscript plus Baseline right-bracket squared equals f Subscript i comma j Superscript 2 Baseline h squared comma EndLayout
Equation (2.5)
StartLayout 1st Row with Label left-parenthesis 2.5 right-parenthesis EndLabel left-bracket left-parenthesis x minus a 1 right-parenthesis Superscript plus Baseline right-bracket squared plus left-bracket left-parenthesis x minus a 2 right-parenthesis Superscript plus Baseline right-bracket squared plus midline-horizontal-ellipsis plus left-bracket left-parenthesis x minus a Subscript n Baseline right-parenthesis Superscript plus Baseline right-bracket squared equals f Subscript i comma j Superscript 2 Baseline h squared EndLayout
Equation (2.6)
StartLayout 1st Row with Label left-parenthesis 2.6 right-parenthesis EndLabel left-parenthesis x minus a 1 right-parenthesis squared plus left-parenthesis x minus a 2 right-parenthesis squared plus midline-horizontal-ellipsis plus left-parenthesis x minus a Subscript p Baseline right-parenthesis squared equals f Subscript i comma j Superscript 2 Baseline h squared and a Subscript p Baseline less-than x overbar less-than-or-equal-to a Subscript p plus 1 Baseline semicolon EndLayout
Equation (2.7)
StartLayout 1st Row with Label left-parenthesis 2.7 right-parenthesis EndLabel StartLayout 1st Row 1st Column left-bracket left-parenthesis u Subscript i comma j Superscript h Baseline 2nd Column minus u Subscript i minus 1 comma j Superscript h Baseline right-parenthesis Superscript plus Baseline right-bracket Superscript 2 Baseline plus left-bracket left-parenthesis u Subscript i comma j Superscript h Baseline minus u Subscript i plus 1 comma j Superscript h Baseline right-parenthesis Superscript plus Baseline right-bracket squared 2nd Row 1st Column Blank 2nd Column plus left-bracket left-parenthesis u Subscript i comma j Superscript h Baseline minus u Subscript i comma j minus 1 Superscript h Baseline right-parenthesis Superscript plus Baseline right-bracket squared plus left-bracket left-parenthesis u Subscript i comma j Superscript h Baseline minus u Subscript i comma j plus 1 Superscript h Baseline right-parenthesis Superscript plus Baseline right-bracket squared equals f Subscript i comma j Superscript 2 Baseline h squared comma 3rd Row 1st Column Blank 2nd Column i equals 2 comma period period period comma upper I minus 1 comma j equals 2 comma period period period comma upper J minus 1 comma EndLayout EndLayout
Equation (2.8)
StartLayout 1st Row with Label left-parenthesis 2.8 right-parenthesis EndLabel StartLayout 1st Row left-bracket left-parenthesis u Subscript i comma j Superscript h Baseline minus u Subscript x max Superscript h Baseline right-parenthesis Superscript minus Baseline right-bracket squared plus left-bracket left-parenthesis u Subscript i comma j Superscript h Baseline minus u Subscript y max Superscript h Baseline right-parenthesis Superscript minus Baseline right-bracket squared equals f Subscript i comma j Superscript 2 Baseline h squared comma 2nd Row i equals 2 comma period period period comma upper I minus 1 comma j equals 2 comma period period period comma upper J minus 1 comma EndLayout EndLayout
Lemma 3.1

Let x overbar be the solution to equation Equation2.5. We have

StartLayout 1st Row with Label left-parenthesis 3.1 right-parenthesis EndLabel 1 greater-than-or-equal-to StartFraction partial-differential x overbar Over partial-differential a 1 EndFraction greater-than-or-equal-to StartFraction partial-differential x overbar Over partial-differential a 2 EndFraction greater-than-or-equal-to midline-horizontal-ellipsis greater-than-or-equal-to StartFraction partial-differential x overbar Over partial-differential a Subscript n Baseline EndFraction greater-than-or-equal-to 0 EndLayout

and

StartLayout 1st Row with Label left-parenthesis 3.2 right-parenthesis EndLabel StartFraction partial-differential x overbar Over partial-differential a 1 EndFraction plus StartFraction partial-differential x overbar Over partial-differential a 2 EndFraction plus midline-horizontal-ellipsis plus StartFraction partial-differential x overbar Over partial-differential a Subscript n Baseline EndFraction equals 1 period EndLayout
Proposition 3.2

In the Gauss-Seidel iteration for the fast sweeping method, the maximum change of u Superscript h at any grid point is less than or equal to the maximum change of u Superscript h at its neighboring points.

Lemma 3.3

The fast sweeping algorithm is monotone in the initial data.

Corollary 3.5

Let u Superscript left-parenthesis k right-parenthesis and v Superscript left-parenthesis k right-parenthesis be two numerical solutions at the k -th iteration of the fast sweeping algorithm. Let double-vertical-bar dot double-vertical-bar Subscript normal infinity be the maximum norm. We have

(1) double-vertical-bar u Superscript left-parenthesis k right-parenthesis Baseline minus v Superscript left-parenthesis k right-parenthesis Baseline double-vertical-bar Subscript normal infinity Baseline less-than-or-equal-to double-vertical-bar u Superscript left-parenthesis k minus 1 right-parenthesis Baseline minus v Superscript left-parenthesis k minus 1 right-parenthesis Baseline double-vertical-bar Subscript normal infinity ,

(2) 0 less-than-or-equal-to max Underscript i comma j Endscripts left-brace u Subscript i comma j Superscript left-parenthesis k right-parenthesis Baseline minus u Subscript i comma j Superscript left-parenthesis k plus 1 right-parenthesis Baseline right-brace less-than-or-equal-to max Underscript i comma j Endscripts left-brace u Subscript i comma j Superscript left-parenthesis k minus 1 right-parenthesis Baseline minus u Subscript i comma j Superscript left-parenthesis k right-parenthesis Baseline right-brace .

Theorem 3.6

The iterative solution by the fast sweeping algorithm converges monotonically to the solution of the discretized system.

Theorem 4.1

For a single data point normal upper Gamma equals StartSet bold x 0 EndSet , the numerical solution, u Superscript h Baseline left-parenthesis bold x comma bold x 0 right-parenthesis , of the fast sweeping method converges in 2 Superscript n sweeps in upper R Superscript n and satisfies

d left-parenthesis bold x comma bold x 0 right-parenthesis less-than-or-equal-to u Superscript h Baseline left-parenthesis bold x comma bold x 0 right-parenthesis less-than-or-equal-to d left-parenthesis bold x comma bold x 0 right-parenthesis plus upper O left-parenthesis StartAbsoluteValue h log h EndAbsoluteValue right-parenthesis comma

where d left-parenthesis bold x comma bold x 0 right-parenthesis is the distance function to bold x 0 .

Equation (4.2)
StartLayout 1st Row with Label left-parenthesis 4.2 right-parenthesis EndLabel StartLayout 1st Row d Subscript x Baseline left-parenthesis x comma y right-parenthesis equals StartFraction x Over StartRoot x squared plus y squared EndRoot EndFraction greater-than 0 comma 0 less-than d Subscript x x Baseline left-parenthesis x comma y right-parenthesis equals StartFraction y squared Over left-parenthesis x squared plus y squared right-parenthesis Superscript 3 slash 2 Baseline EndFraction less-than-or-equal-to StartFraction 1 Over StartRoot x squared plus y squared EndRoot EndFraction comma 2nd Row d Subscript y Baseline left-parenthesis x comma y right-parenthesis equals StartFraction y Over StartRoot x squared plus y squared EndRoot EndFraction greater-than 0 comma 0 less-than d Subscript y y Baseline left-parenthesis x comma y right-parenthesis equals StartFraction x squared Over left-parenthesis x squared plus y squared right-parenthesis Superscript 3 slash 2 Baseline EndFraction less-than-or-equal-to StartFraction 1 Over StartRoot x squared plus y squared EndRoot EndFraction period EndLayout EndLayout
Theorem 4.3

For an arbitrary set of discrete points normal upper Gamma equals StartSet bold x Subscript m Baseline EndSet Subscript m equals 1 Superscript upper M , the numerical solution u Superscript h Baseline left-parenthesis bold x comma normal upper Gamma right-parenthesis by the fast sweeping method after 2 Superscript n sweeps, satisfies

ModifyingAbove u With bar Superscript h Baseline left-parenthesis bold x comma normal upper Gamma right-parenthesis less-than-or-equal-to u Superscript h Baseline left-parenthesis bold x comma normal upper Gamma right-parenthesis less-than-or-equal-to d left-parenthesis bold x comma normal upper Gamma right-parenthesis plus upper O left-parenthesis StartAbsoluteValue h log h EndAbsoluteValue right-parenthesis period
Equation (4.7)
StartLayout 1st Row with Label left-parenthesis 4.7 right-parenthesis EndLabel StartAbsoluteValue d left-parenthesis bold x comma normal upper Gamma overbar right-parenthesis minus d left-parenthesis bold x comma normal upper Gamma right-parenthesis EndAbsoluteValue equals upper O left-parenthesis h right-parenthesis comma EndLayout
Equation (4.8)
StartLayout 1st Row with Label left-parenthesis 4.8 right-parenthesis EndLabel u Superscript h Baseline left-parenthesis bold x comma normal upper Gamma right-parenthesis minus u Superscript h Baseline left-parenthesis bold x comma normal upper Gamma overbar right-parenthesis equals upper O left-parenthesis h right-parenthesis comma for-all bold x comma EndLayout
Example 1

This example shows the interaction of two data points on a simple grid in two dimensions. Five iterations are needed for the fast sweeping algorithm to converge. However, changes after four iterations are upper O left-parenthesis h right-parenthesis no matter what size the grid is as we have tested. In Figure 6.1 we show the results after each iteration on a 7 times 7