LINEBURG


<< . .

 3
( 23)



. . >>

On the Cartesian grid this kink is slightly smeared out, and the derivative
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

<< . .

 3
( 23)



. . >>

Copyright Design by: Sunlight webdesign