LINEBURG

ful (harmless) manner. More important, if the failure of an equation that is

will have a ¬nite value. In fact, consideration of the possible placement of

true in a general sense causes overall degradation of the numerical method,

sample points shows that the value of the derivative lies in the interval

then many times a special-case approach can ¬x the problem. For example,

[’1, 1]. Thus, nothing special needs to be done for kinks. In the worst-case

when calculating the normals using equation (1.2) in the last chapter, we

scenario, the gradient vanishes at a kink, and remedies for this were already

treated the special case where the denominator | φ| was identically zero by

addressed in the last chapter.

randomly choosing the normal direction. The numerical methods outlined

In two spatial dimension we replace the implicit function φ(x) = x2 +

in Part II of this book are based on vanishing viscosity solutions that guar-

y 2 ’ 1 with the signed distance function φ(x) = x2 + y 2 ’ 1 in order to

antee reasonable behavior even at the occasional kink where a derivative

implicitly represent the unit circle ‚„¦ = {x | |x| = 1}. Here | φ| = 1 for

fails to exist.

all x = 0, and a multidimensional kink exists at the single point x = 0.

Again, on our Cartesian grid the kink will be rounded out slightly and will

not pose a problem. However, this numerical smearing of the kink makes

2.4 Examples | φ| = 1 locally. That is, locally φ is no longer a signed distance function,

and one has to take care when applying formulas that assume | φ| = 1.

In the last chapter we used φ(x) = x2 ’ 1 as an implicit representation Luckily, this does not generally lead to catastrophic di¬culties. In fact,

of ‚„¦ = {’1, 1}. A signed distance function representation of these same these kinks mostly exist away from the zero isocontour, which is the region

points is φ(x) = |x|’1, as shown in Figure 2.2. The signed distance function of real interest in interface calculations.

2.5. Geometry and Calculus Toolboxes 21 22 2. Signed Distance Functions

In three spatial dimensions we replace the implicit function φ(x) = x2 + which should not be confused with x, which is the size of a Cartesian

y 2 + z 2 ’ 1 with the signed distance function φ(x) = x2 + y 2 + z 2 ’ 1 grid cell. While this overuse of notation may seem confusing at ¬rst, it is

very common and usually clari¬ed from the context in which it is used.

in order to represent the surface of the unit sphere ‚„¦ = {x | |x| = 1}

Note the simplicity of equation (2.7) as compared to equation (1.8).

implicitly. Again, the multidimensional kink at x = 0 will be smeared out

Obviously, there is much to be gained in simplicity and e¬ciency in us-

on our Cartesian grid.

ing signed distance functions. However, one should be relatively cautious,

In all three examples there was a kink at a single point. This is somewhat

since smeared-out kinks will generally have | φ| = 1, so that equa-

misleading in general. For example, consider the one-dimensional example

tions (2.5) and (2.6) do not accurately de¬ne the normal and the curvature.

φ(x) = |x| ’ 1 again, but in two spatial dimensions, where we write φ(x) =

In fact, when using numerical approximations, one will not generally obtain

|x| ’ 1. Here, the interface consists of the two lines x = ’1 and x = 1,

and the interior region is „¦’ = {x | |x| < 1}. In this example every point | φ| = 1, so equations (2.5) and (2.6) will not generally be accurate. There

are many instances of the normal or the curvature appearing in a set of

along the line x = 0 has a kink in the x direction; i.e., there is an entire

equations when these quantities may not actually be needed or desired. In

line of kinks. Similarly, in three spatial dimensions φ(x) = |x| ’ 1 implicitly

fact, one may actually prefer the gradient of φ (i.e., φ) instead of the nor-

represents the two planes x = ’1 and x = 1. In this case every point on

mal. Similarly, one may prefer the Laplacian of φ (i.e., φ) instead of the

the two-dimensional plane x = 0 has a kink in the x direction; i.e., there is

curvature. In this sense one should always keep equations (2.5) and (2.6)

an entire plane of kinks. All of these kinks will be numerically smeared out

in mind when performing numerical calculations. Even if they are not gen-

on our Cartesian grid, and we need not worry about the derivative being

erally true, they have the potential to make the calculations more e¬cient

unde¬ned. However, locally | φ| = 1 numerically.

and even better behaved in some situations.

The multidimensional delta function in equation (1.19) can be rewritten

as

2.5 Geometry and Calculus Toolboxes

ˆ

δ(x) = δ(φ(x)) (2.8)

Boolean operations for signed distance functions are similar to those for

using the one-dimensional delta function δ(φ). The surface integral in

general implicit functions. If φ1 and φ2 are two di¬erent signed distance

equation (1.21) then becomes

functions, then φ(x) = min(φ1 (x), φ2 (x)) is the signed distance func-

tion representing the union of the interior regions. The function φ(x) =

f (x)δ(φ(x))dx, (2.9)

max(φ1 (x), φ2 (x)) is the signed distance function, representing the intersec- „¦

tion of the interior regions. The complement of the set de¬ned by φ1 (x) has

where the | φ| term has been omitted.

signed distance function φ(x) = ’φ1 (x). Also, φ(x) = max(φ1 (x), ’φ2 (x))

is the signed distance function for the region de¬ned by subtracting the

interior of φ2 from the interior of φ1 .

As mentioned in the last chapter, we would like our implicit function to be

as smooth as possible. It turns out that signed distance functions, especially

those where the kinks have been numerically smeared, are probably the

best candidates for implicit representation of interfaces. This is because

| φ| = 1 everywhere except near the smoothed-out kinks. This simpli¬es

many of the formulas from the last chapter by removing the normalization

constants. Equation (1.2) simpli¬es to

N= φ (2.5)

for the local unit normal. Equation (1.7) simpli¬es to

κ= φ (2.6)

for the curvature, where φ is the Laplacian of φ de¬ned as

φ = φxx + φyy + φzz , (2.7)

Part II 3

Level Set Methods Motion in an Externally Generated

Velocity Field

Level set methods add dynamics to implicit surfaces. The key idea that

started the level set fanfare was the Hamilton-Jacobi approach to numer-

ical solutions of a time-dependent equation for a moving implicit surface.

This was ¬rst done in the seminal work of Osher and Sethian [126]. In the

following chapters we will discuss this seminal work along with many of the

auxiliary equations that were developed along the way, including a general

numerical approach for Hamilton-Jacobi equations.

In the ¬rst chapter we discuss the basic convection equation, otherwise

known as the “level set equation.” This moves an implicit surface in an ex-

3.1 Convection

ternally generated velocity ¬eld. In the following chapter we discuss motion

by mean curvature, emphasizing the parabolic nature of this equation as op-

Suppose that the velocity of each point on the implicit surface is given as

posed to the underlying hyperbolic nature of the level set equation. Then,

V (x); i.e., assume that V (x) is known for every point x with φ(x) = 0.

in the following chapter we introduce the general concept of Hamilton-

Jacobi equations, noting that basic convection is a simple instance of this Given this velocity ¬eld V = u, v, w , we wish to move all the points on

general framework. In the next chapter we discuss the concept of a sur- the surface with this velocity. The simplest way to do this is to solve the

face moving normal to itself. The next two chapters address two of the ordinary di¬erential equation (ODE)

core level set equations and give details for obtaining numerical solutions

dx

in the Hamilton-Jacobi framework. Speci¬cally, we discuss reinitialization = V (x) (3.1)

dt

to a signed distance function and extrapolation of a quantity away from

or across an interface. After this, we discuss a recently developed particle

for every point x on the front, i.e., for all x with φ(x) = 0. This is the

level set method that hybridizes the Eulerian level set approach with La-

Lagrangian formulation of the interface evolution equation. Since there are

grangian particle-tracking technology. Finally, we wrap up this part of the

generally an in¬nite number of points on the front (except, of course, in one

book with a brief discussion of codimension-two (and higher) objects.

spatial dimension), this means discretizing the front into a ¬nite number

of pieces. For example, one could use segments in two spatial dimensions

or triangles in three spatial dimensions and move the endpoints of these

segments or triangles. This is not so hard to accomplish if the connectiv-

ity does not change and the surface elements are not distorted too much.

Unfortunately, even the most trivial velocity ¬elds can cause large distor-

tion of boundary elements (segments or triangles), and the accuracy of the

method can deteriorate quickly if one does not periodically modify the dis-

cretization in order to account for these deformations by smoothing and

regularizing inaccurate surface elements. The interested reader is referred to

[174] for a rather recent least-squares-based smoothing scheme for damping

26 3. Motion in an Externally Generated Velocity Field 3.1. Convection 27

mesh-instabilities due to deforming elements. Examples are given in both On a Cartesian grid it can be slightly complicated to implement equa-

two and three spatial dimensions. Reference [174] also discusses the use of tion (3.2) if the velocity ¬eld is de¬ned only on the interface itself. So,

a mesh-re¬nement procedure to maintain some degree of regularity as the as with the embedding of φ, we usually write equation (3.2) using the

interface deforms. Again, without these special procedures for maintaining assumption that the velocity ¬eld V is not only de¬ned on the interface

both smoothness and regularity, the interface can deteriorate to the point where φ(x) = 0, but is de¬ned o¬ the interface as well. Often V will be

where numerical results are so inaccurate as to be unusable. In addition naturally de¬ned on the entire computational domain „¦, but for numerical

to dealing with element deformations, one must decide how to modify the purposes it is usually su¬cient to have V de¬ned on a band containing

interface discretization when the topology changes. These surgical meth- the interface. The bandwidth varies based on the numerical method used

ods of detaching and reattaching boundary elements can quickly become to obtain approximate solutions to equation (3.2). When V is already de-

rather complicated. Reference [174] outlines some of the details involved in ¬ned throughout all of „¦ nothing special need be done. However, there

a single “surgical cut” of a three-dimensional surface. The use of the La- are interesting examples where V is known only on the interface, and one

grangian formulation of the interface motion given in equation (3.1) along must extend its de¬nition to (at least) a band about the interface in order

with numerical techniques for smoothing, regularization, and surgery are to solve equation (3.2). We will discuss the extension of a velocity o¬ the

collectively referred to as front tracking methods. A seminal work in the interface in Chapter 8.

¬eld of three-dimensional front tracking is [168], and the interested reader Embedding V on our Cartesian grid introduces the same sampling issues

is referred to [165] for a current state-of-the-art review. that we faced in Chapter 1 when we embedded the interface “ as the zero

In order to avoid problems with instabilities, deformation of surface level set of the function φ. For example, suppose we were given a velocity

elements, and complicated surgical procedures for topological repair of in- ¬eld V that is identically zero in all of „¦ except on the interface, where

terfaces, we use our implicit function φ both to represent the interface V = 1, 0, 0 . Then the exact solution is an interface moving to the right

and to evolve the interface. In order to de¬ne the evolution of our implicit with speed 1. However, since most (if not all) of the Cartesian grid points

function φ we use the simple convection (or advection) equation will not lie on the interface, most of the points on our Cartesian mesh have

V identically equal to zero, causing the V · φ term in equation (3.2) to

φt + V · φ = 0, (3.2)

vanish. This in turn implies that φt = 0 almost everywhere, so that the

where the t subscript denotes a temporal partial derivative in the time interface mostly (or completely if no points happen to fall on the interface)

variable t. Recall that is the gradient operator, so that incorrectly sits still. This di¬cult issue can be recti¬ed in part by placing

some conditions on the velocity ¬eld V . For example, if we require that

V· φ = uφx + vφy + wφz .

V be continuous near the interface, then this rules out our degenerate

This partial di¬erential equation (PDE) de¬nes the motion of the interface example.

where φ(x) = 0. It is an Eulerian formulation of the interface evolution, Restricting V to the set of continuous functions generally does not alle-

since the interface is captured by the implicit function φ as opposed to being viate our sampling problems. Suppose, for example, that the above velocity

tracked by interface elements as was done in the Lagrangian formulation. ¬eld was equal to 1, 0, 0 on the interface, zero outside a band of thickness

Equation (3.2) is sometimes referred to as the level set equation; it was > 0 surrounding the interface, and smooth in between. We can choose V

introduced for numerical interface evolution by Osher and Sethian [126]. It as smooth as we like by de¬ning it appropriately in the band of thickness

is also a quite popular equation in the combustion community, where it is surrounding the interface. The di¬culty arises when is small compared to

known as the G-equation given by x. If is small enough, then almost every grid point will lie outside the

band where V = 0. Once again, we will (mostly) compute an interface that

Gt + V · G = 0, (3.3)

incorrectly sits still. In fact, even if is comparable to x, the numerical

where the G(x) = 0 isocontour is used to represent implicitly the reac- solution will have signi¬cant errors. In order to resolve the velocity ¬eld, it

tion surface of an evolving ¬‚ame front. The G-equation was introduced by is necessary to have a number of grid points within the thickness band

Markstein [110], and it is used in the asymptotic analysis of ¬‚ame fronts in surrounding the interface. That is, we need x to be signi¬cantly smaller

instances where the front is thin enough to be considered a discontinuity. than the velocity variation (which scales like ) in order get a good ap-

The interested reader is referred to Williams [173] as well. Lately, numeri- proximation of the velocity near the interface. Since x needs to be much

cal practitioners in the combustion community have started using level set smaller than , we desire a relatively large to minimize the variation in

methods to ¬nd numerical solutions of equation (3.3) in (obviously) the the velocity ¬eld.

same manner as equation (3.2).

28 3. Motion in an Externally Generated Velocity Field 3.2. Upwind Di¬erencing 29

3.2 Upwind Di¬erencing

Given a velocity ¬eld V and the notion (discussed above) that minimizing

its variation is good for treating the sampling problem, there is an obvious

choice of V that gives both the correct interface motion and the least vari- Once φ and V are de¬ned at every grid point (or at least su¬ciently close

ation. First, since the values of V given on the interface dictate the correct to the interface) on our Cartesian grid, we can apply numerical methods

interface motion, these cannot be changed, regardless of the variation. In to evolve φ forward in time moving the interface across the grid. At some

some sense, the spatial variation of the velocity on the interface dictates point in time, say time tn , let φn = φ(tn ) represent the current values

how many Cartesian grid points will be needed to accurately predict the in- of φ. Updating φ in time consists of ¬nding new values of φ at every grid

terface motion. If we cannot resolve the tangential variation of the interface point after some time increment t. We denote these new values of φ by

velocity with our Cartesian grid, then it is unlikely that we can calculate φn+1 = φ(tn+1 ), where tn+1 = tn + t.

a good approximation to the interface motion. Second, the velocity o¬ the A rather simple ¬rst-order accurate method for the time discretization

interface has nothing to do with the correct interface motion. This is true of equation (3.2) is the forward Euler method given by

even if the velocity o¬ the interface is inherited from some underlying phys-

φn+1 ’ φn

+Vn· φn = 0,

ical calculation. Only the velocity of the interface itself contains any real (3.4)

t

information about the interface propagation. Otherwise, one would have

no hope of using the Lagrangian formulation, equation (3.1), to calculate where V n is the given external velocity ¬eld at time tn , and φn evaluates

the interface motion. In summary, the velocity variation tangential to the the gradient operator using the values of φ at time tn . Naively, one might

interface dictates the interface motion, while the velocity variation normal evaluate the spatial derivatives of φ in a straightforward manner using equa-

to the interface is meaningless. Therefore, the minimum variation in the tion (1.3), (1.4), or (1.5). Unfortunately, this straightforward approach will

velocity ¬eld can be obtained by restricting the interface velocity V to be fail. One generally needs to exercise great care when numerically discretiz-

constant in the direction normal to the interface. This generally makes the ing partial di¬erential equations. We begin by writing equation (3.4) in

velocity multivalued, since lines normal to the interface will eventually in- expanded form as

tersect somewhere away from the interface (if the interface has a nonzero

φn+1 ’ φn

+ un φ n + v n φ n + w n φ n = 0

curvature). Alternatively, the velocity V (x) at a point x can be set equal to (3.5)

x y z

t

the interface velocity V (xC ) at the interface point xC closest to the point x.

and address the evaluation of the un φn term ¬rst. The techniques used to

While this doesn™t change the value of the velocity on the interface, it makes x

approximate this term can then be applied independently to the v n φn and

the velocity o¬ the interface approximately constant in the normal direction y

wn φn terms in a dimension-by-dimension manner.

local to the interface. In Chapter 8 we will discuss numerical techniques for z

For simplicity, consider the one-dimensional version of equation (3.5),

constructing a velocity ¬eld de¬ned in this manner.

De¬ning the velocity V equal to the interface velocity at the closest in- φn+1 ’ φn

+ un φn = 0, (3.6)

terface point xC is a rather ingenious idea. In the appendix of [175], Zhao x

t

et al. showed that a signed distance function tends to stay a signed dis-

where the sign of un indicates whether the values of φ are moving to the

tance function if this closest interface point velocity is used to advect the

right or to the left. Since un can be spatially varying, we focus on a speci¬c

interface. A number of researchers have been using this specially de¬ned ve-

grid point xi , where we write

locity ¬eld because it usually gives superior results over velocity ¬elds with

needlessly more spatial variation. Chen, Merriman, Osher, and Smereka n+1

φi ’ φni

+ un (φx )n = 0, (3.7)

[43] published the ¬rst numerical results based on this specially designed i i

t

velocity ¬eld. The interested reader is referred to the rather interesting

where (φx )i denotes the spatial derivative of φ at the point xi . If ui > 0,

work of Adalsteinsson and Sethian [1] as well.

the values of φ are moving from left to right, and the method of charac-

The velocity ¬eld given in equation (3.2) can come from a number of ex-

teristics tells us to look to the left of xi to determine what value of φ will

ternal sources. For example, when the φ(x) = 0 isocontour represents the

land on the point xi at the end of a time step. Similarly, if ui < 0, the

interface between two di¬erent ¬‚uids, the interface velocity is calculated

values of φ are moving from right to left, and the method of characteristics

using the two-phase Navier-Stokes equations. This illustrates that the in-

implies that we should look to the right to determine an appropriate value

terface velocity more generally depends on both space and time and should

of φi at time tn+1 . Clearly, D’ φ (from equation (1.4)) should be used to

be written as V (x, t), but we occasionally omit the x dependence and more

approximate φx when ui > 0. In contrast, D+ φ cannot possibly give a good

often the t dependence for brevity.

30 3. Motion in an Externally Generated Velocity Field 3.3. Hamilton-Jacobi ENO 31

approximation, since it fails to contain the information to the left of xi that although

dictates the new value of φi . Similar reasoning indicates that D + φ should

max{|V |}

be used to approximate φx when ui < 0. This method of choosing an ap- t =± (3.11)

min{ x, y, z}

proximation to the spatial derivatives based on the sign of u is known as

upwind di¬erencing or upwinding. Generally, upwind methods approximate

is also quite popular. More details on consistency, stability, and conver-

derivatives by biasing the ¬nite di¬erence stencil in the direction where the

gence can be found in basic textbooks on the numerical solution of partial

characteristic information is coming from.

di¬erential equations; see, for example, [157].

We summarize the upwind discretization as follows. At each grid point,

Instead of upwinding, the spatial derivatives in equation (3.2) could be

de¬ne φ’ as D’ φ and φ+ as D+ φ. If ui > 0, approximate φx with φ’ . If

x x x approximated with the more accurate central di¬erencing. Unfortunately,

ui < 0, approximate φx with φ+ . When ui = 0, the ui (φx )i term vanishes,

x simple central di¬erencing is unstable with forward Euler time discretiza-

and φx does not need to be approximated. This is a ¬rst-order accurate

tion and the usual CFL conditions with t ∼ x. Stability can be achieved

discretization of the spatial operator, since D ’ φ and D+ φ are ¬rst-order

by using a much more restrictive CFL condition with t ∼ ( x)2 , al-

accurate approximations of the derivative; i.e., the errors are O( x).

though this is too computationally costly. Stability can also be achieved

The combination of the forward Euler time discretization with the up-

by using a di¬erent temporal discretization, e.g., the third-order accurate

wind di¬erence scheme is a consistent ¬nite di¬erence approximation to

Runge-Kutta method (discussed below). A third way of achieving stabil-

the partial di¬erential equation (3.2), since the approximation error con-

ity consists in adding some arti¬cial dissipation to the right-hand side of

verges to zero as t ’ 0 and x ’ 0. According to the Lax-Richtmyer

equation (3.2) to obtain

equivalence theorem a ¬nite di¬erence approximation to a linear partial

di¬erential equation is convergent, i.e., the correct solution is obtained as φt + V · φ = µ φ, (3.12)

t ’ 0 and x ’ 0, if and only if it is both consistent and stable. Stability

where the viscosity coe¬cient µ is chosen proportional to x, i.e., µ ∼ x,

guarantees that small errors in the approximation are not ampli¬ed as the

so that the arti¬cial viscosity vanishes as x ’ 0, enforcing consistency

solution is marched forward in time.

for this method. While all three of these approaches stabilize central di¬er-

Stability can be enforced using the Courant-Friedreichs-Lewy condition

encing, we instead prefer to use upwind methods, which draw on the highly

(CFL condition), which asserts that the numerical waves should propagate

sucessful technology developed for the numerical solution of conservation

at least as fast as the physical waves. This means that the numerical wave

laws.

speed of x/ t must be at least as fast as the physical wave speed |u|,

i.e., x/ t > |u|. This leads us to the CFL time step restriction of

3.3 Hamilton-Jacobi ENO

x

t< , (3.8)

max{|u|}

The ¬rst-order accurate upwind scheme described in the last section can

where max{|u|} is chosen to be the largest value of |u| over the entire be improved upon by using a more accurate approximation for φ’ and φ+ .

x x

Cartesian grid. In reality, we only need to choose the largest value of |u| on The velocity u is still used to decide whether φ’ or φ+ is used, but the

x x

the interface. Of course, these two values are the same if the velocity ¬eld is approximations for φ’ or φ+ can be improved signi¬cantly.

x x

de¬ned as the velocity of the closest point on the interface. Equation (3.8) In [81], Harten et al. introduced the idea of essentially nonoscillatory

is usually enforced by choosing a CFL number ± with (ENO) polynomial interpolation of data for the numerical solution of con-

servation laws. Their basic idea was to compute numerical ¬‚ux functions

max{|u|}

using the smoothest possible polynomial interpolants. The actual numerical

t =± (3.9)

x implementation of this idea was improved considerably by Shu and Osher

in [150] and [151], where the numerical ¬‚ux functions were constructed

and 0 < ± < 1. A common near-optimal choice is ± = 0.9, and a common

directly from a divided di¬erence table of the pointwise data. In [126],

conservative choice is ± = 0.5. A multidimensional CFL condition can be

Osher and Sethian realized that Hamilton-Jacobi equations in one spatial

written as

dimension are integrals of conservation laws. They used this fact to extend

the ENO method for the numerical discretization of conservation laws to

|u| |v| |w|

t max + + = ±, (3.10) Hamilton-Jacobi equations such as equation (3.2). This Hamilton-Jacobi

x y z

32 3. Motion in an Externally Generated Velocity Field 3.4. Hamilton-Jacobi WENO 33

of φ+ . In other words, ¬rst-order accurate polynomial interpolation is

ENO (HJ ENO) method allows one to extend ¬rst-order accurate upwind x

di¬erencing to higher-order spatial accuracy by providing better numerical exactly ¬rst-order upwinding. Improvements are obtained by including

approximations to φ’ or φ+ . the Q2 (xi ) and Q3 (xi ) terms in equation (3.18), leading to second- and

x x

Proceeding along the lines of [150] and [151], we use the smoothest pos- third-order accuracy, respectively.

1

sible polynomial interpolation to ¬nd φ and then di¬erentiate to get φx . As Looking at the divided di¬erence table and noting that Dk+1/2 φ was

is standard with Newton polynomial interpolation (see any undergraduate chosen for ¬rst-order accuracy, we have two choices for the second-order

numerical analysis text, e.g., [82]), the zeroth divided di¬erences of φ are accurate correction. We could include the next point to the left and use

2 2

de¬ned at the grid nodes and de¬ned by Dk φ, or we could include the next point to the right and use Dk+1 φ. The

key observation is that smooth slowly varying data tend to produce small

0

Di φ = φ i (3.13)

numbers in divided di¬erence tables, while discontinuous or quickly vary-

ing data tend to produce large numbers in divided di¬erence tables. This

at each grid node i (located at xi ). The ¬rst divided di¬erences of φ are

is obvious in the sense that the di¬erences measure variation in the data.

de¬ned midway between grid nodes as

2 2

Comparing |Dk φ| to |Dk+1 φ| indicates which of the polynomial interpolants

0 0

Di+1 φ ’ Di φ

1

has more variation. We would like to avoid interpolating near large varia-

Di+1/2 φ = , (3.14)

x tions such as discontinuities or steep gradients, since they cause overshoots

where we are assuming that the mesh spacing is uniformly x. Note that in the interpolating function, leading to numerical errors in the approxi-

Di’1/2 φ = (D’ φ)i and Di+1/2 φ = (D+ φ)i , i.e., the ¬rst divided dif-

1 1 2 2 2

mation of the derivative. Thus, if |Dk φ| ¤ |Dk+1 φ|, we set c = Dk φ and

2

ferences, are the backward and forward di¬erence approximations to the k = k ’ 1; otherwise, we set c = Dk+1 φ and k = k. Then we de¬ne

derivatives. The second divided di¬erences are de¬ned at the grid nodes as

Q2 (x) = c(x ’ xk )(x ’ xk+1 ), (3.21)

1 1

’

Di+1/2 φ Di’1/2 φ

so that

2

Di φ = , (3.15)

2x

Q2 (xi ) = c (2(i ’ k) ’ 1) x (3.22)

while the third divided di¬erences

is the second-order accurate correction to the approximation of φx in

D 2 φ ’ Di φ

2

= i+1 equation (3.18). If we stop here, i.e., omitting the Q3 term, we have a

3

Di+1/2 φ (3.16)

3x second-order accurate method for approximating φ’ and φ+ . Note that

x x

k has not yet been used. It is de¬ned below for use in calculating the

are de¬ned midway between the grid nodes.

third-order accurate correction.

The divided di¬erences are used to reconstruct a polynomial of the form

Copyright Design by: Sunlight webdesign