Fast Sweeping Algorithms for a Class of
Hamilton-Jacobi Equations
Yen-Hsi Richard Tsai
∗†
, Li-Tien Cheng
‡§
, Stanley Osher
¶†
, and Hong-Kai Zhao
k
Abstract
We derive a Godunov-type numerical flux for the class of strictly convex,
homogeneous Hamiltonians that includes H(p, q) =
p
ap
2
+ bq
2
− 2cpq,
c
2
< ab. We combine our Godunov numerical fluxes with simple Gauss-
Seidel type iterations for solving the corresponding Hamilton-Jacobi Equa-
tions. The resulting algorithm is fast since it does not require a sorting strat-
egy as found, e.g., in the fast marching method. In addition, it provides a way
to compute solutions to a class of HJ equations for which the conventional
fast marching method is not applicable. Our experiments indicate conver-
gence after a few iterations, even in rather difficult cases.
1 Introduction
Hamilton-Jacobi equations have a rich pool of applications, ranging from those of
optimal control theory, geometrical optics, to essentially any problem that needs the
(weighted) distance function [13]. Examples include crystal growth, ray tracing,
etching, robotic motion planning , and computer vision. Solutions of these types
of equations usually develop singularities in their derivatives, and thus, the unique
viscosity solution [5] is sought.
∗
Department of Mathematics and PACM, Princeton University, Princeton, New Jersey 08644.
Email: ytsai@math.princeton.edu
†
Research supported by ONR N00014-97-1-0027, DARPA/NSF VIP grant NSF DMS 9615854
and ARO DAAG 55-98-1-0323
‡
Department of Mathematics, UCSD, La Jolla, CA 92093-0112. Email: lcheng@math.ucsd.edu
§
Research supported by NSF Grant #0112413 and #0208449
¶
Department of Mathematics, UCLA, Los Angeles, California 90095. Email: sjo@math.ucla.edu
k
Department of Mathematics, UCI, Irvine, CA 92697-3875. Email: zhao@math.ucla.edu
1
In this article, we focus on the class of time independent Hamilton-Jacobi Equa-
tions with Dirichlet boundary condition:
H(x, ∇u) = r(x), u|
Γ
= 0;
H(x, p) are strictly convex, non-negative, and lim
λ→0
H(x, λp )=0. We explain
our method using this following important model equation:
H(φ
x
, φ
y
) =
q
aφ
2
x
+ bφ
2
y
− 2cφ
x
φ
y
= r, (1)
where φ : R
2
7→ R continuous and a, b, c and r can either be constants or scalar
functions, in which case H depends also on x, defined on R
2
, satisfying ab > c
2
,
a, b, r > 0. With a = b = 1, and c = 0, we have the standard eikonal equation
for which many numerical methods have been developed. This equation has the
essential features of HJ equations with convex Hamiltonians so that we can easily
explain our algorithm, and is general enough that fast marching is not applicable.
In the following subsections, we will review some of the solution methods for the
eikonal equation since it forms the motivation of our work. We then present a fast
Gauss-Seidel type iteration method for equation (1) which utilizes a monotone up-
wind Godunov flux for the Hamiltonian. We show numerically that this algorithm
can be applied directly to equations of the above type with variable coefficients.
1.1 Solving Eikonal Equations
In geometrical optics [9], the eikonal equation
q
φ
2
x
+ φ
2
y
= r(x, y) (2)
is derived from the leading term in an asymptotic expansion
e
iω(φ(x,y)−t)
∞
X
j=0
A
j
(x, y, t)(iω)
−j
of the wave equation:
w
tt
− c
2
(x, y)(w
xx
+ w
yy
) = 0,
where r(x, y) = 1/|c(x, y)|, is the function of slowness. The level sets of the
solution φ can thus be interpreted as the first arrival time of the wave front that is
initially Γ. It can also be interpreted as the “distance” function to Γ.
2
We first restrict our attention for now to the case in which r = 1. Let Γ be a closed
subset of R
2
. It can be shown easily that the distance function defined by
d(x) = dist (x, Γ) := min
p∈Γ
||x − p||, x = (x, y) ∈ R
2
,
is the viscosity solution to equation (2) with the boundary condition
φ(x, y) = 0 for (x, y) ∈ Γ.
Rouy and Tourin [19] proved the convergence to the viscosity solution of an iter-
ative method solving equation (2) with the Godunov Hamiltonian approximating
||∇φ||. The Godunov Hamiltonian function can be written in the following form:
H
G
(p
−
, p
+
, q
−
, q
+
) =
q
max{p
+
−
, p
−
+
}
2
+ max{q
+
−
, q
−
+
}
2
(3)
where p
±
= D
x
±
φ
i,j
, q
±
= D
y
±
φ
i,j
, D
x
±
φ
i,j
= ±(φ
i±1,j
−φ
i,j
)/h and accordingly
for D
y
±
φ
i,j
, and x
+
= max(x, 0), x
−
= − min(x, 0).
Osher [12] provided a link to time dependent eikonal equations by proving that the
t-level set of φ(x, y) is the zero level set of the viscosity solution of the evolution
equation at time t
ψ
t
= ||∇ψ||
with appropriate initial conditions. In fact, the same is true for a very general
class of Hamilton-Jacobi equations (see [12]). As a consequence, one can try to
solve the time-dependent equation by the level set formulation [16] with high order
approximations on the partial derivatives [8][17]. Crandall and Lions proved that
the discrete solution obtained with a consistent, monotone Hamiltonian converges
to the desired viscosity solution [4].
Tsitsiklis [24] combined heap sort with a variant of the classical Dijkstra algorithm
to solve the steady state equation of the more general problem
||∇φ|| = r(x).
This was later rederived in [22] and also reported in [7]. It has become known as the
fast marching method whose complexity is O(N log(N)), where N is the number
of grid points. Osher and Helmsen [14] have extended the fast marching type
method to somewhat more general Hamilton-Jacobi equations. We will comment
on this in a following section.
3
1.2 Anisotropic Eikonal Equation
We return to the Hamiltonian in question: H(p, q) =
p
ap
2
+ bq
2
− 2cpq. Writing
the quadratic form as
ap
2
+ bq
2
− 2cpq =
p q
a −c
−c b
p
q
,
it is easy to see that we can diagonalize the symmetric matrix in the middle of the
equation for our previously noted choices of a, b, c and find a coordinate system
ξ-η such that after rescaling, the Hamiltonian becomes
H(˜p, ˜q) =
p
˜p
2
+ ˜q
2
.
The eigensystem of the above matrix defines the anisotropy. Indeed, the authors in
[14] proposed to solve the constant coefficient equation (1) by first transforming it
to equation (2) in the ξ-η coordinate system.
This anisotropy occurs in fields such as ray tracing in special media, e.g. crystals,
in which there are “preferred” directions. Furthermore, we will see that it can be a
result of considering the geodesic distance function on a manifold M that is defined
as the graph of a smooth function f.
Let φ be the distance function such that
φ(x, y) = min
γ⊂M
Z
γ
ds,
and γ connects the point (x, y) with the set Γ ⊂ M. The minimizing curve is called
the geodesic, and φ the distance function to Γ on M. Moreover, φ solves
||P
∇ψ
∇φ||
2
= 1, φ|
Γ
= 0, (4)
where ψ(x, y, z) = f(x, y) − z, and the projection operator [3]
P
∇ψ
= I −
∇ψ
N
∇ψ
||∇ψ||
2
,
which projects a vector onto a plane whose normal is parallel to ∇ψ. Using the
fact that P
∇ψ
is a projection operator, a simple calculation shows that
||P
∇ψ
∇φ||
2
= (1−
f
2
x
f
2
x
+ f
2
y
+ 1
)φ
2
x
+(1−
f
2
y
f
2
x
+ f
2
y
+ 1
)φ
2
y
−2
f
x
f
y
f
2
x
+ f
2
y
+ 1
φ
x
φ
y
.
(5)
4
This is clearly of the form of Hamiltonians that we are interested in. We will apply
our algorithm to compute the geodesic distance later in this paper.
There are other approaches that are designed to compute distances on manifolds.
[10] provided an algorithm to compute the geodesic distance on triangulated mani-
folds. Barth [2] uses the Discontinuous Galerkin Method to find distance on graphs
of functions that are represented by spline functions. In [3], the authors embed
the manifold as the zero level set of a Lipschitz continuous function and solve
the corresponding time-dependent eikonal equation (4) in the embedding space.
As we have mentioned in the previous subsection, the zero level set of the time-
dependent eikonal equation at time t
1
is the t
1
-level set of the solution to the sta-
tionary eikonal equations (see [12]). In [11], the authors adopted the standard fast
marching method to solve the isotropic eikonal equation in a thin band of thickness
, that encloses the manifold M, and proved that the restriction of the solution to
M converges to the geodesic distance as goes to 0. In [20, 21], the authors pro-
vide an ordered upwind method to solve a general class of static Hamilton-Jacobi
equations. We will comment on their method in a later subsection.
1.3 Osher’s Fast Marching Criteria
Since the fast marching method is by now well known, we will not give much detail
on its implementation in this paper. In general, this involves a sorting procedure
and the solution of
H
G
(D
x
−
φ
i,j
, D
x
+
φ
i,j
, D
y
−
φ
i,j
, D
y
+
φ
i,j
) = 1 (6)
for φ
ij
in terms of its four neighboring values. More precisely, the heap sort strat-
egy of the fast marching method requires a monotone update sequence. The up-
dated value of a grid node has to be greater than or equal to those of the grid nodes
used to form the finite difference stencil. This amounts to the condition
pH
p
+ qH
q
≥ 0,
which dictates that the solution is non-decreasing along the characteristics. How-
ever, if we use one sided upwind finite difference approximations for partial deriva-
tives of φ on a Cartesian grid, it is equivalent to demanding that the partial deriva-
tives of φ (i.e. p and q) and their corresponding components of the characteristics
directions (i.e. dx/dt and dy/dt) have the same sign. Since dx/dt = H
p
and
dy/dt = H
q
, we have the stricter Osher’s fast marching criterion:
pH
p
≥ 0, qH
q
≥ 0. (7)
5
![Figure 1: The 1-level set of H and the box I[p-,p+]xI[q-,q+]](/figures/figure-1-the-1-level-set-of-h-and-the-box-i-p-p-xi-q-q-1k2axtq9.webp)
