LINEBURG


<< . .

 5
( 23)



. . >>

where φ is the temperature and b is the thermal conductivity. The heat
stepping with the CFL condition in equation (4.7) as well.
equation is the most basic equation of the parabolic model.
The stringent O(( x)2 ) time-step restriction resulting from the forward
When φ is a signed distance function, bκ| φ| and b φ are identical, and
Euler time discretization can be alleviated by using an ODE solver with
either of these can be used to calculate the right-hand side of equation (4.5).
a larger stability region, e.g., an implicit method. For example, ¬rst-order
However, once this right-hand side is combined with a forward Euler time
4.3. Convection-Di¬usion Equations 45 46 4. Motion Involving Mean Curvature

accurate backward Euler time stepping applied to equation (4.6) gives and the two can be used interchangeably if one maintains a signed distance
approximation for φ o¬ the interface. These equations can be solved using
φn+1 ’ φn the upwind methods from the last chapter on the V · φ term and cen-
= b φn+1 , (4.8)
t tral di¬erencing on the parabolic b φ or bκ| φ| term. A TVD RK time
discretization can be used with a time-step restriction of
which has no time step stability restriction on the size of t. This means
that t can be chosen for accuracy reasons alone, and one typically sets
|u| |v| |w| 2b 2b 2b
t = O( x). Note that setting t = O( x) as opposed to t = O(( x)2 ) t + + + + + <1 (4.12)
2 2 2
x y z ( x) ( y) ( z)
lowers the overall accuracy to O( x). This can be improved upon using
the trapezoidal rule satis¬ed at every grid point.
Suppose the O(1) size b term is replaced with an O( x) size term
φn+1 ’ φn φn + φn+1
=b , (4.9) that vanishes as the mesh is re¬ned with x ’ 0. Then equation (4.10)
t 2
becomes
which is O(( t)2 ) in time and thus O(( x)2 ) overall even when t =
φt + V · φ= φ, (4.13)
O( x). This combination of the trapezoidal rule with central di¬erencing
which asymptotically approaches equation (3.2) as ’ 0. The addition of
of a parabolic spatial operator is generally referred to as the Crank-Nicolson
an arti¬cial φ term to the right-hand side of equation (3.2) is called the
scheme.
arti¬cial viscosity method. Arti¬cial viscosity is used by many authors to
The price we pay for the larger time step achieved using either equa-
stabilize a central di¬erencing approximation to the convective φ term in
tion (4.8) or equation (4.9) is that a linear system of equations must
be solved at each time step to obtain φn+1 . Luckily, this is not di¬cult equation (3.2). This arises in computational ¬‚uid dynamics, where terms
given the simple linear structure of φn+1 . Unfortunately, an implicit dis- of the form φ are added to the right-hand side of convective equations
to pick out vanishing viscosity solutions valid in the limit as ’ 0. This
cretization of equation (4.5) requires consideration of the more complicated
nonlinear κn+1 | φn+1 | term. vanishing viscosity picks out the physically correct weak solution when no
classical solution exists, for example in the case of a discontinuous shock
We caution the reader that one cannot substitute equation (4.6) for equa-
tion (4.5) when using an implicit time discretization. Even if φn is initially wave. It is interesting to note that the upwind discretizations discussed
a signed distance function, φn+1 will generally not be a signed distance in the last chapter have numerical truncation errors that serve the same
function after the linear system has been solved. This means that φn+1 purpose as the φ term. First-order accurate upwinding has an intrinsic
is not a good approximation to κn+1 | φn+1 | even though φn may be O( x) arti¬cial viscosity, and the higher-order accurate upwind methods
have intrinsic arti¬cial viscosities with magnitude O(( x)r ), where r is the
exactly equal to κn | φn |. Although we stress (throughout the book) the
order of accuracy of the method.
conceptual simpli¬cations and computational savings that can be obtained
In [146], Sethian suggested an entropy condition that required curves to
when φ is a signed distance function, e.g., replacing N with φ, κ with
¬‚ow into corners, and he provided numerical evidence to show that this
φ, etc., we caution the reader that there is a signi¬cant and important
entropy condition produced the correct weak solution for self-interesting
di¬erence between the two in the case where φ is not a signed distance
curves. Sethian™s entropy condition indicates that κ| φ| is a better form
function.
for the vanishing viscosity than φ for dealing with the evolution of lower-
dimensional interfaces. This concept was rigorized by Osher and Sethian in
[126], where they pointed out that
4.3 Convection-Di¬usion Equations
φt + V · φ = κ| φ| (4.14)
The convection-di¬usion equation is a more natural choice than equation (4.13) for dealing with level set
methods, although these two equations are interchangeable when φ is a
φt + V · φ=b φ (4.10)
signed distance function.
includes both the e¬ects of an external velocity ¬eld and a di¬usive term.
The level set version of this is
φt + V · φ = bκ| φ|, (4.11)
48 5. Hamilton-Jacobi Equations

5.2 Connection with Conservation Laws
Consider the one-dimensional scalar conservation law

5 ut + F (u)x = 0; (5.3)

where u is the conserved quantity and F (u) is the ¬‚ux function. A well-
Hamilton-Jacobi Equations known conservation law is the continuity equation

ρt + (ρu)x = 0 (5.4)

for conservation of mass, where ρ is the density of the material. In compu-
tational ¬‚uid dynamics (CFD), the continuity equation is combined with
equations for conservation of momentum and conservation of energy to ob-
tain the compressible Navier-Stokes equations. When viscous e¬ects are
ignored, the Navier-Stokes equations reduce to the compressible inviscid
Euler equations.
The presence of discontinuities in the Euler equations forces one to con-
sider weak solutions where the derivatives of solution variables, e.g., ρx ,
can fail to exist. Examples include linear contact discontinuities and non-
linear shock waves. The nonlinear nature of shock waves allows them to
develop as the solution progresses forward in time even if the data are ini-
5.1 Introduction tially smooth. The Euler equations may not always have unique solutions,
and an entropy condition is used to pick out the physically correct solu-
In this chapter we discuss numerical methods for the solution of general tion. This is the vanishing viscosity solution discussed in the last chapter.
Hamilton-Jacobi equations of the form For example, vanishing viscosity admits physically consistent rarefaction
waves, ruling out physically inadmissible expansion shocks.
Burgers™ equation
φt + H( φ) = 0, (5.1)
u2
ut + =0 (5.5)
where H can be a function of both space and time. In three spatial 2 x
dimensions, we can write
is a scalar conservation law that possesses many of the interesting non-
linear properties contained in the more complex Euler equations. Burgers™
φt + H(φx , φy , φz ) = 0 (5.2)
equation develops discontinuous shock waves from smooth initial data and
exhibits nonphysical expansion shocks if the vanishing viscosity solution is
as an expanded version of equation (5.1). Convection in an externally gen- not used to force these to become smooth rarefaction waves. Many of the
erated velocity ¬eld (equation (3.2)) is an example of a Hamilton-Jacobi numerical methods developed to solve Burgers™ equation can be extended to
equation where H( φ) = V · φ. The level set equation (equation (4.4)) treat both the one-dimensional and the multidimensional Euler equations
is another example of a Hamilton-Jacobi equation with H( φ) = Vn | φ|. of gas dynamics.
Here Vn can depend on x, t, or even φ/| φ|. Consider the one-dimensional Hamilton-Jacobi equation
The equation for motion by mean curvature (equation (4.5)) is not a
φt + H(φx ) = 0, (5.6)
Hamilton-Jacobi-type equation, since the front speed depends on the sec-
ond derivatives of φ. Hamilton-Jacobi equations depend on (at most) the
which becomes
¬rst derivatives of φ, and these equations are hyperbolic. The equation for
motion by mean curvature depends on the second derivatives of φ and is (φx )t + H(φx )x = 0 (5.7)
parabolic.
5.3. Numerical Discretization 49 50 5. Hamilton-Jacobi Equations

after one takes a spatial derivative of the entire equation. Setting u = φx or HJ WENO schemes. For brevity, we discuss the two-dimensional numer-
in equation 5.7 results in ical approximation to H(φx , φy ), noting that the extension to three spatial
dimensions is straightforward. An important class of schemes is that of
ut + H(u)x = 0, (5.8)
monotone schemes. A scheme is monotone when φn+1 as de¬ned in equa-
which is a scalar conservation law; see equation (5.3). Thus, in one spatial tion (5.9) is a nondecreasing function of all the φn . Crandall and Lions
dimension we can draw a direct correspondence between Hamilton-Jacobi proved that these schemes converge to the correct solution, although they
equations and conservation laws. The solution u to a conservation law is are only ¬rst-order accurate. The numerical Hamiltonians associated with
the derivative of a solution φ to a Hamilton-Jacobi equation. Conversely, monotone schemes are important, and examples will be given below.
the solution φ to a Hamilton-Jacobi equation is the integral of a solution u The forward Euler time discretization (equation (5.9)) can be extended to
to a conservation law. This allows us to point out a number of useful facts. higher-order TVD Runge Kutta in a straightforward manner, as discussed
For example, since the integral of a discontinuity is a kink, or discontinuity in Chapter (3). The CFL condition for equation 5.9 is
in the ¬rst derivative, solutions to Hamilton-Jacobi equations can develop
|H1 | |H2 | |H3 |
kinks in the solution even if the data are initially smooth. In addition, t max + + < 1, (5.10)
x y z
solutions to Hamilton-Jacobi equations cannot generally develop a discon-
tinuity unless the corresponding conservation law develops a delta function. where H1 , H2 , and H3 are the partial derivatives of H with respect to φx ,
Thus, solutions φ to equation (5.2) are typically continuous. Furthermore, φy , and φz , respectively. For example, in equation (3.2), where H( φ) =
since conservation laws can have nonunique solutions, entropy conditions V · φ, the partial derivatives of H are H1 = u, H2 = v, and H3 = w. In
are needed to pick out “physically” relevant solutions to equation (5.2) as this case equation (5.10) reduces to equation (3.10). As another example,
well. consider the level set equation (4.4) with H( φ) = Vn | φ|. Here the partial
Viscosity solutions for Hamilton-Jacobi equations were ¬rst proposed by derivatives are slightly more complicated, with H1 = VN φx /| φ|, H2 =
Crandall and Lions [52]. Monotone ¬rst-order accurate numerical meth- VN φy /| φ|, and H3 = VN φz /| φ|, assuming that VN does not depend on
ods were ¬rst presented by Crandall and Lions [53] as well. Later, Osher φx , φy or φz . Otherwise, the partial derivatives can be substantially more
and Sethian [126] used the connection between conservation laws and complicated.
Hamilton-Jacobi equations to construct higher-order accurate “artifact-
free” numerical methods. Even though the analogy between conservation
5.3.1 Lax-Friedrichs Schemes
laws and Hamilton-Jacobi equations fails in multiple spatial dimensions,
many Hamilton-Jacobi equations can be discretized in a dimension by di- ˆ
The ¬rst approximation to H that we consider is the Lax-Friedrichs (LF)
mension fashion. This culminated in [127], where Osher and Shu proposed scheme from [53] given by
a general framework for the numerical solution of Hamilton-Jacobi equa-
φ’ + φ+ φ’ + φ+ φ+ ’ φ’
φ+ ’ φ’
tions using successful methods from the theory of conservation laws. We y y y y
ˆ x x x x
’±x ’±y
H=H , , (5.11)
2 2 2 2
follow [127] below.
where ±x and ±y are dissipation coe¬cients that control the amount of
numerical viscosity. These dissipation coe¬cients
5.3 Numerical Discretization ±x = max |H1 (φx , φy )|, ±y = max |H2 (φx , φy )| (5.12)
are chosen based on the partial derivatives of H.
A forward Euler time discretization of a Hamilton-Jacobi equation can be
The choice of the dissipation coe¬cients in equation (5.12) can be rather
written as
subtle. In the traditional implementation of the LF scheme, the maximum is
φn+1 ’ φn ˆ
+ H n (φ’ , φ+ ; φ’ , φ+ ; φ’ , φ+ ) = 0, (5.9) chosen over the entire computational domain. First, the maximum and min-
x x y y z z
t
imum values of φx are identi¬ed by considering all the values of φ’ and φ+
x x
ˆxxyyzz
where H(φ’ , φ+ ; φ’ , φ+ ; φ’ , φ+ ) is a numerical approximation of H(φx , φy , x min max
on the Cartesian mesh. Then one can identify the interval I = [φx , φx ].
ˆ A similar procedure is used to de¬ne I y = [φy , φmax ]. The coe¬cients
min
φz ). The function H is called a numerical Hamiltonian, and it is required y
ˆ ±x and ±y are set to the maximum possible values of |H1 (φx , φy )| and
to be consistent in the sense that H(φx , φx ; φy , φy ; φz , φz ) = H(φx , φy , φz ).
Recall that spatial derivatives such as φ’ are discretized with either ¬rst- |H2 (φx , φy )|, respectively, with φx ∈ I x and φy ∈ I y . Although it is oc-
x
casionally di¬cult to evaluate the maximum values of |H1 | and |H2 |, it is
order accurate one-sided di¬erencing or the higher-order accurate HJ ENO
5.3. Numerical Discretization 51 52 5. Hamilton-Jacobi Equations

φ± and φ± at xi,j using the HJ WENO scheme. This type of scheme has
straightforward to do so in many instances. For example, in equation (3.2), x y
both H1 = u and H2 = v are independent of φx and φy , so ±x and ±y can been referred to as a Stencil Lax-Friedrichs (SLF) scheme, since it deter-
be set to the maximum values of |u| and |v| on the Cartesian mesh. mines the dissipation coe¬cient using only the neighboring grid points that
Consider evaluating ±x and ±y for equation (4.4) where H1 = VN φx /| φ| are part of the stencil used to determine φx and φy . An alternative to the
and H2 = VN φy /| φ|, recalling that these are the partial derivatives only if dimension-by-dimension neighborhoods is to use the 49 grid points in the
VN is independent φx and φy . It is somewhat more complicated to evaluate rectangle with diagonal corners at xi’3,j’3 and xi+3,j+3 to determine ±.
±x and ±y in this case. When φ is a signed distance function with | φ| = 1 This idea of searching only locally to determine the dissipation coe¬-
(or ≈ 1 numerically), we can simplify to H1 = VN φx and H2 = VN φy . cients can be taken a step further. The Local Lax-Friedrichs (LLF) scheme
These functions can still be somewhat tricky to work with if VN is spatially proposed for conservation laws by Shu and Osher [151] does not look at
varying. But in the special case that VN is spatially constant, the maximum any neighboring grid points when calculating the dissipation coe¬cients in
values of |H1 | and |H2 | can be determined by considering only the endpoints a given direction. In [127], Osher and Shu interpreted this to mean that
±x is determined at each grid point using only the values of φ’ and φ+
of Ix and Iy , respectively. This is true because H1 and H2 are monotone x x
at that speci¬c grid point to determine the interval I x . The interval I y is
functions of φx and φy , respectively. In fact, when VN is spatially constant,
H1 = VN φx /| φ| and H2 = VN φy /| φ| are straightforward to work with as still determined in the LF or SLF manner (in the SLF case we rename LLF
as SLLF). Similarly, ±y uses an interval I y , de¬ned using only the values
well. The function H1 achieves a maximum when |φx | is as large as possible
of φ’ and φ+ at the grid point in question while I x is still determined in
and |φy | is as small as possible. Thus, only the endpoints of I x and I y need y y
be considered; note that we use φy = 0 when the endpoints of I y di¬er the LF or SLF fashion. Osher and Shu [127] also proposed the Local Lo-
in sign. Similar reasoning can be used to ¬nd the maximum value of |H2 |. cal Lax-Friedrichs (LLLF) scheme with even less numerical dissipation. At
each grid point I x is determined using the values of φ’ and φ+ at that grid
One way to treat a spatially varying VN is to make some estimates. For x x
point; I y is determined using the values of φ’ and φ+ at that grid point;
example, since |φx |/| φ| ¤ 1 for all φx and φy , we can bound |H1 | ¤ |VN |. y y
A similar bound of |H2 | ¤ |VN | holds for |H2 |. Thus, both ±x and ±y can and then these intervals are used to determine both ± and ±y . When H is
x

separable, i.e., H(φx , φy ) = H x (φx ) + H y (φy ), LLLF reduces to LLF, since
be set to the maximum value of |VN | on the Cartesian mesh. The price
±x is independent of φy , and ±y is independent of φx . When H is not sep-
we pay for using bounds to choose ± larger than it should be is increased
numerical dissipation. That is, while the numerical method will be stable arable, LLF and LLLF are truly distinct schemes. In practice, LLF seems
and give an accurate solution as the mesh is re¬ned, some details of this to work better than any of the other options. LF and SLF are usually too
solution may be smeared out and lost on a coarser mesh. dissipative, while LLLF is usually not dissipative enough to overcome the
Since increasing ± increases the amount of arti¬cial dissipation, decreas- problems introduced by using the centrally averaged approximation to φx
ing the quality of the solution, it is bene¬cial to chose ± as small as possible and φy in evaluating H in equation (5.11). Note that LLF is a monotone
without inducing oscillations or other nonphysical phenomena into the so- scheme.
ˆ
lution. In approximating Hi,j at a grid point xi,j on a Cartesian mesh,
it then makes little sense to do a global search to de¬ne the intervals I x
and I y . In particular, consider the simple convection equation (3.2) where 5.3.2 The Roe-Fix Scheme
±x = max |u| and ±y = max |v|. Suppose that some region had relatively
As discussed above, choosing the appropriate amount of arti¬cial dissipa-
small values of |u| and |v|, while another region had relatively large values.
tion to add to the centrally evaluated H in equation (5.11) can be tricky.
Since the LF method chooses ±x as the largest value of |u| and ±y as the
Therefore, it is often desirable to use upwind-based methods with built-in
largest value of |v|, the same values of ± will be used in the region where
arti¬cial dissipation. For conservation laws, Shu and Osher [151] proposed
the velocities are small as is used in the region where the velocities are
using Roe™s upwind method along with an LLF entropy correction at sonic
large. In the region where the velocities are large, the large values of ± are
points where entropy-violating expansion shocks might form. The added
required to obtain a good solution. But in the region where the velocities
dissipation from the LLF entropy correction forces the expansion shocks
are small, these large values of ± produce too much numerical dissipation,
to develop into continuous rarefaction waves. The method was dubbed Roe
wiping out small features of the solution. Thus, it is advantageous to use
Fix (RF) and it can be written for Hamilton-Jacobi equations (see [127])
only the grid points su¬ciently close to xi,j in determining ±. A rule of
as
thumb is to include the grid points from xi’3,j to xi+3,j in the x-direction
φ+ ’ φ’
φ+ ’ φ’
and from xi,j’3 to xi,j+3 in the y-direction in the local search neighborhood y y
ˆ x x
H = H φx , φy ’ ± x y
’± , (5.13)
for determining ±. This includes all the grid nodes that are used to evaluate 2 2
5.3. Numerical Discretization 53 54 5. Hamilton-Jacobi Equations

where ±x and ±y are usually set identically to zero in order to remove formulas and use these cheaper approximations to decide whether or not
the numerical dissipation terms. In the RF scheme, I x and I y are initially upwinding or LLF will be used. After making this decision, the higher-
determined using only the nodal values for φ± and φ± as in the LLLF order accurate HJ WENO (or HJ ENO) method can be used to compute
x y
the necessary values of φ± and φ± used in the numerical discretization,
scheme. In order to estimate the potential for upwinding, we look at the x y
partial derivatives H1 and H2 . If H1 (φx , φy ) has the same sign (either obtaining the usual high-order accuracy. Sonic points rarely occur in prac-
always positive or always negative) for all φx ∈ I x and all φy ∈ I y , we tice, and this strategy reduces the use of the costly HJ WENO method by
know which way information is ¬‚owing and can apply upwinding. Similarly, approximately a factor of two.
if H2 (φx , φy ) has the same sign for all φx ∈ I x and φy ∈ I y , we can
upwind this term as well. If both H1 and H2 do not change sign, we upwind
5.3.3 Godunov™s Scheme
completely, setting both ±x and ±y to zero. If H1 > 0, information is
¬‚owing from left to right, and we set φx = φ’ . Otherwise, H1 < 0, and we In [74], Godunov proposed a numerical method that gives the exact solu-
x
set φx = φ+ . Similarly, H2 > 0 indicates φy = φ’ , and H2 < 0 indicates tion to the Riemann problem for one-dimensional conservation laws with
x y
φy = φ+ . piecewise constant initial data. The multidimensional Hamilton-Jacobi
y
If either H1 or H2 changes sign, we are in the vicinity of a sonic point formulation of this scheme can be written as
where the eigenvalue (in this case H1 or H2 ) is identically zero. This signi¬es
ˆ
H = extx exty H (φx , φy ) , (5.14)
a potential di¬culty with nonunique solutions, and arti¬cial dissipation
is needed to pick out the physically correct vanishing viscosity solution. as was pointed out by Bardi and Osher [12]. This is the canonical monotone
We switch from the RF scheme to the LLF scheme to obtain the needed scheme. De¬ning our intervals I x and I y in the LLLF manner using only
arti¬cial viscosity. If there is a sonic point in only one direction, i.e., x or y, the values of φ± and φ± at the grid node under consideration, we de¬ne
x y
it makes little sense to add damping in both directions. Therefore, we look extx and exty as follows. If φ’ < φ+ , then extx H takes on the minimum
x x
for sonic points in each direction and add damping only to the directions value of H for all φx ∈ I x . If φ’ > φ+ , then extx H takes on the maximum
x x
that have sonic points. This is done using the I x and I y de¬ned as in the value of H for all φx ∈ I x . Otherwise, if φ’ = φ+ , then extx H simply plugs
x x
LLF method. That is, we switch from the LLLF de¬ned intervals used φ’ (= φ+ ) into H for φx . Similarly, if φ’ < φ+ , then exty H takes on the
x x y y
above to initially look for sonic points to the larger LLF intervals that are minimum value of H for all φy ∈ I y . If φ’ > φ+ , then exty H takes on the
y y
even more likely to have sonic points. We proceed as follows. If H1 (φx , φy ) maximum value of H for all φy ∈ I y . Otherwise, if φ’ = φ+ , then exty H
y y
y
x
does not change sign for all φx ∈ ILLF and all φy ∈ ILLF , we set φx equal to simply plugs φ’ (= φ+ ) into H for φy . In general, extx exty H = exty extx H,
y y
either φ’ or φ+ depending on the sign of H1 . In addition, we set ±x to zero so di¬erent versions of Godunov™s method are obtained depending on the
x x
to remove the arti¬cial dissipation in the x-direction. At the same time, this order of operations. However, in many cases, including when H is separable,
means that a sonic point must have occurred in H2 , so we use an LLF-type extx exty H = exty extx H so this is not an issue.
method for the y-direction, setting φy = (φ’ + φ+ )/2 and choosing ±y as Although Godunov™s method can sometimes be di¬cult to implement,
y y
dictated by the LLF scheme. A similar algorithm is executed if H2 (φx , φy ) there are times when it is straightforward. Consider equation (3.2) for mo-
y
x
does not change sign for all φx ∈ ILLF and φy ∈ ILLF . Then φy is set to tion in an externally generated velocity ¬eld. Here, we can consider the
either φ’ or φ+ , depending on the sign of H2 ; ±y is set to zero; and an LLF x and y directions independently, since H is separable with extx exty H =
y y
method is used in the x-direction, setting φx = (φ’ + φ+ )/2 while choosing extx (uφx ) + exty (vφy ). If φ’ < φ+ , we want the minimum value of uφx .
x x x x
±x as dictated by the LLF scheme. If both H1 and H2 change sign, we have Thus, if u > 0, we use φ’ , and if u < 0, we use φ+ . If u = 0, we obtain
x x
sonic points in both directions and proceed with the standard LLF scheme uφx = 0 regardless of the choice of φx . On the other hand, if φ’ > φ+ ,
x x
at that grid point. we want the maximum value of uφx . Thus, if u > 0, we use φ’ , and ifx
With the RF scheme, upwinding in the x-direction dictates that either u < 0, we use φ+ . Again, u = 0 gives uφx = 0. Finally, if φ’ = φ+ , then
x x x

φx or φ+ be used, but not both. Similarly, upwinding in the y-direction uses uφx is uniquely determined. This can be summarized as follows. If u > 0,
x
either φ’ or φ+ , but not both. Since evaluating φ± and φ± using high-order use φ’ ; if u < 0, use φ+ ; and if u = 0, set uφx = 0. This is identical to
y y x y x x
accurate HJ ENO or HJ WENO schemes is rather costly, it seems wasteful the standard upwind di¬erencing method described in Chapter 3. That is,
to do twice as much work in these instances. Unfortunately, one cannot for motion in an externally generated velocity ¬eld, Godunov™s method is
determine whether upwinding can be used (as opposed to LLF) without identical to simple upwind di¬erencing.
computing φ± and φ± . In order to minimize CPU time, one can compute
x y
φ± and φ± using the ¬rst-order accurate forward and backward di¬erence
x y
56 6. Motion in the Normal Direction

1 1 1

0.5 0.5 0.5

0 0 0

6 ’0.5 ’0.5 ’0.5

’1 ’1 ’1
’1 0 1 ’1 0 1 ’1 0 1
Motion in the Normal Direction 1 1 1

0.5 0.5 0.5

0 0 0

’0.5 ’0.5 ’0.5

’1 ’1 ’1
’1 0 1 ’1 0 1 ’1 0 1

1 1 1

0.5 0.5 0.5

0 0 0

’0.5 ’0.5 ’0.5

’1 ’1 ’1
’1 0 1 ’1 0 1 ’1 0 1
6.1 The Basic Equation
Figure 6.1. Evolution of a star-shaped interface as it moves normal to itself in
In this chapter we discuss the motion of an interface under an internally the outward direction.
generated velocity ¬eld for constant motion in the normal direction. This
velocity ¬eld is de¬ned by V = aN or Vn = a, where a is a constant. The
corresponding level set equation (i.e., equation (4.4)) is Thus, if φn is initially a signed distance function (with | φn | = 1), it stays
a distance function (with | φ| = 1) for all time.
φt + a| φ| = 0, (6.1)
When the initial data constitute a signed distance function, forward Euler
time stepping reduces to solving the ordinary di¬erential equation φt = ’a
where a can be of either sign. When a > 0 the interface moves in the
independently at every grid point. Since a is a constant, this forward Euler
normal direction, and when a < 0 the interface moves opposite the normal
time stepping gives the exact solution up to round-o¬ error (i.e., there is no
direction. When a = 0 this equation reduces to the trivial φt = 0, where
truncation error). For example, consider a point where φ = φo > 0, which is

<< . .

 5
( 23)



. . >>

Copyright Design by: Sunlight webdesign