Eikonal equation explained
An eikonal equation (from Greek εἰκών, image[1] [2]) is a non-linear first-order partial differential equation that is encountered in problems of wave propagation.
The classical eikonal equation in geometric optics is a differential equation of the form
where
lies in an open subset of
,
is a positive function,
denotes the
gradient, and
is the Euclidean norm. The function
is given and one seeks solutions
.In the context of
geometric optics, the function
is the
refractive index of the medium.
More generally, an eikonal equation is an equation of the form
where
is a function of
variables. Here the function
is given, and
is the solution.If
, then equation becomes .
Eikonal equations naturally arise in the WKB method[3] and the study of Maxwell's equations. Eikonal equations provide a link between physical (wave) optics and geometric (ray) optics.
One fast computational algorithm to approximate the solution to the eikonal equation is the fast marching method.
History
The term "eikonal" was first used in the context of geometric optics by Heinrich Bruns.[4] However, the actual equation appears earlier in the seminal work of William Rowan Hamilton on geometric optics.[5]
Physical interpretation
Continuous shortest-path problems
Suppose that
is an open set with suitably smooth boundary
. The solution to the eikonal equation
\left|\nablau(x)\right|=
forx\in\Omega\subsetRn,
u(x)=q(x)forx\in\partial\Omega
can be interpreted as the minimal amount of time required to travel from
to
, where
f:\bar{\Omega}\to(0,+infty)
is the speed of travel, and
q:\partial\Omega\to[0,+infty)
is an exit-time penalty. (Alternatively this can be posed as a minimal cost-to-exit by making the right-side
and
an exit-cost penalty.)
In the special case when
, the solution gives the
signed distance from
.
[6] By assuming that
exists at all points, it is easy to prove that
corresponds to a time-optimal control problem using Bellman's optimality principle and a Taylor expansion.
[7] Unfortunately, it is not guaranteed that
exists at all points, and more advanced techniques are necessary to prove this. This led to the development of
viscosity solutions in the 1980s by
Pierre-Louis Lions and
Michael G. Crandall,
[8] and Lions won a
Fields Medal for his contributions.
Electromagnetic potential
The physical meaning of the eikonal equation is related to the formula
where
is the electric field strength, and
is the electric potential. There is a similar equation for velocity potential in fluid flow and temperature in heat transfer. The physical meaning of this equation in the electromagnetic example is that any charge in the region is pushed to move at right angles to the lines of constant potential, and along lines of force determined by the field of the
E vector and the sign of the charge.
Ray optics and electromagnetism are related by the fact that the eikonal equation gives a second electromagnetic formula of the same form as the potential equation above where the line of constant potential has been replaced by a line of constant phase, and the force lines have been replaced by normal vectors coming out of the constant phase line at right angles. The magnitude of these normal vectors is given by the square root of the relative permittivity. The line of constant phase can be considered the edge of one of the advancing light waves (wavefront). The normal vectors are the rays the light is traveling down in ray optics.
Computational algorithms
Several fast and efficient algorithms to solve the eikonal equation have been developed since the 1990s. Many of these algorithms take advantage of algorithms developed much earlier for shortest path problems on graphs with nonnegative edge lengths.[9] These algorithms take advantage of the causality provided by the physical interpretation and typically discretize the domain using a mesh[10] [11] [12] [13] or regular grid[14] [15] and calculate the solution at each discretized point.Eikonal solvers on triangulated surfaces were introduced in.
Sethian's fast marching method (FMM) was the first "fast and efficient" algorithm created to solve the Eikonal equation. The original description discretizes the domain
into a regular grid and "marches" the solution from "known" values to the undiscovered regions, precisely mirroring the logic of
Dijkstra's algorithm. If
is discretized and has
meshpoints, then the computational complexity is
where the
term comes from the use of a heap (typically binary). A number of modifications can be prescribed to FMM since it is classified as a label-setting method. In addition, FMM has been generalized to operate on general meshes that discretize the domain.
Label-correcting methods such as the Bellman–Ford algorithm can also be used to solve the discretized Eikonal equation also with numerous modifications allowed (e.g. "Small Labels First" [16] or "Large Labels Last" [17]). Two-queue methods have also been developed[18] that are essentially a version of the Bellman-Ford algorithm except two queues are used with a threshold used to determine which queue a gridpoint should be assigned to based on local information.
Sweeping algorithms such as the fast sweeping method (FSM)[19] are highly efficient for solving Eikonal equations when the corresponding characteristic curves do not change direction very often. These algorithms are label-correcting but do not make use of a queue or heap, and instead prescribe different orderings for the gridpoints to be updated and iterate through these orderings until convergence. Some improvements were introduced such as "locking" gridpoints during a sweep if does not receive an update, but on highly refined grids and higher-dimensional spaces there is still a large overhead due to having to pass through every gridpoint. Parallel methods have been introduced that attempt to decompose the domain and perform sweeping on each decomposed subset. Zhao's parallel implementation decomposes the domain into
-dimensional subsets and then runs an individual FSM on each subset.
[20] Detrixhe's parallel implementation also decomposes the domain, but parallelizes each individual sweep so that processors are responsible for updating gridpoints in an
-dimensional
hyperplane until the entire domain is fully swept.
[21] Hybrid methods have also been introduced that take advantage of FMM's efficiency with FSM's simplicity. For example, the Heap Cell Method (HCM) decomposes the domain into cells and performs FMM on the cell-domain, and each time a "cell" is updated FSM is performed on the local gridpoint-domain that lies within that cell. A parallelized version of HCM has also been developed.[22]
Numerical approximation
For simplicity assume that
is discretized into a uniform grid with spacings
and
in the x and y directions, respectively.
2D approximation on a Cartesian grid
Assume that a gridpoint
has value
. A first-order scheme to approximate the partial derivatives is
U,
U,0\right)2+
U,
U,0\right)2 =
where
ux(xij) ≈
U=
and uy(xij) ≈
U=
.
Due to the consistent, monotone, and causal properties of this discretization it is easy to show that if
and
and
then
\right)2+\left(
\right)2=
which can be solved as a quadratic. In the limiting case of
, this reduces to
Uij=
+
\sqrt{(UX+U
-2
+
-
\right)}.
This solution will always exist as long as
is satisfied and is larger than both,
and
, as long as
.
If
, a lower-dimensional update must be performed by assuming one of the partial derivatives is
:
Uij=min\left(UX+
,UY+
\right).
n-D approximation on a Cartesian grid
Assume that a grid point
has value
. Repeating the same steps as in the
case we can use a first-order scheme to approximate the partial derivatives. Let
be the minimum of the values of the neighbors in the
directions, where
is a
standard unit basis vector. The approximation is then
Solving this quadratic equation for
yields:
If the discriminant in the square root is negative, then a lower-dimensional update must be performed (i.e. one of the partial derivatives is
).
If
then perform the one-dimensional update
If
then perform an
dimensional update using the values
\{U1,\ldots,Un\}\setminus\{Ui\}
for every
and choose the smallest.
Mathematical description
An eikonal equation is one of the form
u(0,x')=u0(x'),forx=(x1,x')
The plane
can be thought of as the initial condition, by thinking of
as
We could also solve the equation on a subset of this plane, or on a curved surface, with obvious modifications.
c2|\nablaxu|2=|\partialtu|2
, where
and
. In geometric optics, the eikonal equation describes the phase fronts of waves. Under reasonable hypothesis on the "initial" data, the eikonal equation admits a local solution, but a global smooth solution (e.g. a solution for all time in the geometrical optics case) is not possible. The reason is that
caustics may develop. In the geometrical optics case, this means that wavefronts cross.
We can solve the eikonal equation using the method of characteristics. One must impose the "non-characteristic" hypothesis
along the initial hypersurface
, where
H =
H(
x,
p) and
p = (
p1,...,
pn) is the variable that gets replaced by ∇
u. Here
x = (
x1,...,
xn) = (
t,
x′).
First, solve the problem
,
. This is done by defining curves (and values of
on those curves) as
=\nabla\xiH(x(s),\xi(s)),
=-\nablaxH(x(s),\xi(s)).
x(0)=x0, \xi(x(0))=\nablau(x(0)).
Note that even before we have a solution
, we know
for
due to our equation for
.That these equations have a solution for some interval
follows from standard ODE theorems (using the non-characteristic hypothesis). These curves fill out an
open set around the plane
. Thus the curves define the value of
in an open set about our initial plane. Once defined as such it is easy to see using the chain rule that
\partialsH(x(s),\xi(s))=0
, and therefore
along these curves.
We want our solution
to satisfy
, or more specifically, for every
,
(\nablau)(x(s))=\xi(x(s)).
Assuming for a minute that this is possible, for any solution
we must have
u(x(s))=\nablau(x(s)) ⋅
=\xi ⋅
,
and therefore
u(x(t))=u(x(0))+
\xi(x(s)) ⋅
ds.
In other words, the solution
will be given in a neighborhood of the initial plane by an explicit equation. However, since the different paths
, starting from different initial points may cross, the solution may become multi-valued, at which point we have developed caustics.We also have (even before showing that
is a solution)
\xi(x(t))=\xi(x(0))-
\nablaxH(x(s),\xi(x(s)))ds.
It remains to show that
, which we have defined in a neighborhood of our initial plane, is the gradient of some function
. This will follow if we show that the vector field
is curl free. Consider the first term in the definition of
. This term,
is curl free as it is the gradient of a function. As for the other term, we note
| \partial2 |
\partialxk\partialxj |
H=
| \partial2 |
\partialxj\partialxk |
H.
The result follows.
Applications
See also
Further reading
- Book: Paris . D. T. . Hurd . F. K. . Basic Electromagnetic Theory . McGraw-Hill . 1969 . 383–385 . 0-07-048470-8 .
- Book: Arnold, V. I. . Vladimir Arnold . Lectures on Partial Differential Equations . Springer . 2004 . 2nd . 2–3 . 3-540-40448-1 .
External links
Notes and References
- The Oxford English Dictionary. 2nd ed. 1989. OED Online. Oxford University Press. 4 April 2000 http://dictionary.oed.com/cgi/entry/00292404
- Book: Evans, L. C. . Partial Differential Equations . AMS Graduate Texts in Mathematics . 19 . 93 .
- Book: Dimassi . Mouez . Sjöstrand . Johannes . Spectral asymptotics in the semi-classical limit . London Math. Society Lecture Notes 268. . Cambridge University Press . 1999 . 0-521-66544-2.
- Book: Bruns, Heinrich . Das Eikonal . S. Hirzel . 1895.
- Hamilton . William Rowan . 1828 . Theory of Systems of Rays . Transactions of the Royal Irish Academy . 15 . 69–174.
- Sakai, Takashi. "On Riemannian manifolds admitting a function whose gradient is of constant norm." Kodai Mathematical Journal 19.1 (1996): 39-51.
- Z. . Clawson . A. . Chacon . A. . Vladimirsky . Causal Domain Restriction for Eikonal Equations . SIAM Journal on Scientific Computing . 36 . 5 . A2478–A2505 . 2014 . 1309.2884 . 10.1137/130936531 . 2014SJSC...36A2478C . 17226196 .
- Book: M. . Bardi . I. . Capuzzo-Dolcetta . Optimal Control and Viscosity Solutions of Hamilton–Jacobi–Bellman Equations . Boston . Birkhäuser . 1997 . 0-8176-3640-4 .
- A. . Chacon . A. . Vladimirsky . Fast two-scale methods for eikonal equations . SIAM Journal on Scientific Computing . 34 . 2 . A547–A578 . 2012 . 10.1137/10080909X . 1110.6220 . 2012SJSC...34A.547C . 6404391 .
- R. . Kimmel . J. A. . Sethian . Computing Geodesic Paths on Manifolds . Proceedings of the National Academy of Sciences . 95 . 15 . 8431–8435 . 1998 . 10.1073/pnas.95.15.8431 . 9671694 . 21092 . 1998PNAS...95.8431K . free .
- A. M. . Bronstein . M. M. . Bronstein . R. . Kimmel . Weighted distance maps computation on parametric three-dimensional manifolds . Journal of Computational Physics . 225 . 1 . 771–784 . 2007 . 10.1016/j.jcp.2007.01.009 . 2007JCoPh.225..771B .
- J. A. . Sethian . A. . Vladimirsky . Fast methods for the Eikonal and related Hamilton–Jacobi equations on unstructured meshes . Proc. Natl. Acad. Sci. USA . 97 . 11 . 5699–5703 . 2000 . 10.1073/pnas.090060097 . 10811874 . 18495 . 2000PNAS...97.5699S . free .
- D. S. . Yershov . S. M. . LaValle . Simplicial Dijkstra and A* Algorithms: From Graphs to Continuous Spaces . Advanced Robotics . 26 . 17 . 2065–2085 . 2012 . 10.1080/01691864.2012.729559 . 17573584 .
- J. A. . Sethian . A Fast Marching Level Set Method for Monotonically Advancing Fronts . Proc. Natl. Acad. Sci. . 93 . 4 . 1591–1595 . 1996 . 10.1073/pnas.93.4.1591 . 11607632 . 39986 . 1996PNAS...93.1591S . free .
- J. N. . Tsitsiklis . Efficient algorithms for globally optimal trajectories . IEEE Trans. Autom. Control . 40 . 9 . 1528–1538 . 1995 . 10.1109/9.412624 . 1721.1/3340 . free .
- D. P. . Bertsekas . A Simple and Fast Label Correcting Algorithm for Shortest Paths . Networks . 23 . 8 . 703–709 . 1993 . 10.1002/net.3230230808 . 1721.1/3256 . free .
- D. P. . Bertsekas . F. . Guerriero . R. . Musmanno . Parallel Asynchronous Label Correcting Methods for Shortest Paths . Journal of Optimization Theory and Applications . 88 . 297–320 . 1996 . 2 . 10.1007/BF02192173 . 13172492 . 1721.1/3390 . free .
- S. . Bak . J. . McLaughlin . D. . Renzi . Some improvements for the fast sweeping method . SIAM Journal on Scientific Computing . 32 . 5 . 2010 . 2853–2874 . 10.1137/090749645 . 2010SJSC...32.2853B .
- H. . Zhao . A fast sweeping method for eikonal equations . Math. Comp. . 74 . 2004 . 250 . 603–627 . 10.1090/S0025-5718-04-01678-3 . free .
- H. . Zhao . Parallel Implementations of the Fast Sweeping Method . J. Comput. Math. . 25 . 4 . 421–429 . 2007 . 43693378 .
- M. . Detrixhe . F. . Gibou . C. . Min . A parallel fast sweeping method for the Eikonal equation . Journal of Computational Physics . 237 . 46–55 . 2013 . 10.1016/j.jcp.2012.11.042 . 2013JCoPh.237...46D .
- A. . Chacon . A. . Vladimirsky . A parallel two-scale method for Eikonal equations . SIAM Journal on Scientific Computing . 37 . 1 . A156–A180 . 2015 . 10.1137/12088197X . 1306.4743 . 2015SJSC...37A.156C .