LINEBURG

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

Copyright Design by: Sunlight webdesign