2.2.

Methods for Obtaining FD Expressions

There are several methods, and we will look at a few:
1) Taylor series expansion the most common, but purely mathematical.
2) Polynomial fitting or interpolation the most general ways. Taylor series is a subset of this method. Interpolation takes us
back to the M.O.C. and thus has a more physical interpretation.
3) Control volume approach also called finite volume (FV) we solve the equations in integral rather than differential form.
Popular in engineering where complex geometries and coordinate transformations are involved. For Cartesian grids, simplest
FV methods FD.
We will look at only the first two approaches.

2.2.1. Taylor Series Expansion Method

Recall the definition of a derivative:
u
x

= lim
x0 , y0

x 0

u ( x0 + x, y0 ) u ( x0 , y0 )
x

The Taylor series approach works backwards want to approximate u/x by a discrete difference, i.e., for finite x.
Given u(x0 , y0 ), we can write a Taylor series expansion for u(x0 +x, y0 ), as
u ( x0 + x) = u ( x0 ) + x

u
x

+
x0 , y0

( x) 2 2u
2
2! x

+
x0 , y0

(x )3 3u
3
3! x

( x) n nu
n
n ! x
n= 0

+ ... =
x0 , y0

This expression is exact if we retain all terms!

The grid mesh or stencil looks like:

Here, ui+1,j = u( x0 +x, y0 )

ui,j+1 = u( x0 , y0 +y)
etc.
If we solve for u/x from the Taylor series, we have
u
x

=
x0 , y0

u ( x0 + x, y0 ) u ( x0 , y0 )
x

x 2 u
2
2! x

x0 , y0

( x) 2 3u
3
3! x

+ ...

(1)

x0 , y0

The first term on the RHS is simply the slope of the function u(x,y,), using the current and the points to the right.

Therefore,

u
x

u ( x0 + x, y0 ) u ( x0 , y0 )
x

x0 , y0

+ O( x ) .

One can also use the value at x0 x instead to get

u
x

=
x0 , y0

u ( x0 , y0 ) u( x0 x , y0 )
x

x 2 u
2
2! x

x0 , y0

( x) 2 3u
3
3! x

+ ...

(2)

x0 , y0

Both (1) and (2) provides an expression for u/x, but numerically the answers will be different.
(1) is called a forward difference
(2) is called a backward difference

Consider a 1-D example:

If we have the equation

u
u
+c
= 0, where c > 0,
t
x

we might want to use

u
u u
+ c i i 1 0,
t
x

which is called Upwind Difference. Alternatively, in

u
u u
+ c i +1 i 0,
t
x

Downstream Difference is used.

Upstream difference is better than downstream difference for this problem, because for this pure advection problem, signals move
from upstream (left) to downstream (right). The value of u at point i at a future time should be influenced by the values of u upstream,
not downstream.

Think of it in terms of characteristics:

No information is coming from the +x direction. This of course depends on the sign of c. If C<0, then upstream means the left side of
the current point.
Note that upstream or upstream-biased schemes are usually better choices in CFD (obviously we are talking about hyperbolic
equations that represent propagations!).
We will see later, that 1st-order downstream difference is absolutely unstable for the above problem.
We can get another discrete approximation to u/x by adding (1) and (2):
u
x

=
x0 , y0

ui+1 ui1 ( x ) 2 3u

3
2 x
3! x

+ ...

(3)

x0 , y0

This is called Centered Difference. Note, it doesn't even use value of u at the current point i. It approximates the slope using two
neighboring points:

Note that this will not be accurate if u has very sharp gradient (). Note also that the extra term in (3) are different they control the
accuracy.
We can build many different approximations to the derivatives through linear combinations of various T.S. expansions.

For high-order derivatives, we follow the same approach.

Consider 2 u/x2 . The Taylor series gives
u ( x + x ) = u ( x ) + x

u (x ) 2 2u ( x ) 3 3u
+
+
+ ...
x
2! x 2
3! x 3

(4)

2u
u
from (4), but we don't want to have the unknown
.
2
x
x
We can replace it with one of the earlier approximations to it, or make use of the following:
We can solve for

u ( x x ) = u ( x ) x

u ( x ) 2 2u ( x )3 3u
+

+ ...
x
2! x 2
3! x3

and add (4) and (5) and solve for

(5)

2u
:
x 2

2u u ( x + x) 2u ( x) + u ( x x)
=
+ O( x 2 )
2
2
x
( x)
2u
which is a centered difference for 2 .
x
Truncation Error
The high-order terms (H.O.T.) of O(xn ) are called the Truncation Error, and are a measure of the error associated with representing a
PDE by a truncated T.S. we can't retain all terms, so we neglect the terms of order O(x2 ) and above, for example.
PDE = FDE +

where is the Truncation Error.

If we change the nature of the H.O.T., by using a different approximation, the accuracy of the F.D. expression will also change.
It is important to understand the impact of , and this is the subject of our computer problem #4.
From the above, it's clear that we want 0 when PDE = FDE. Otherwise, we have a problem (consistency)!
Let x 0, then should 0 => our discrete system approaches continuum and our FDE PDE.
Order of the F.D. Approximation.
The power to which the leading discrete interval in is raised is called the Order of the F.D. Approximation.
Example:

ui ui 1
is first-order accurate in space
x

because =

x 2 u
+ ... = O( x1)
2
2! x

ui+1 2ui + ui1

is second order accurate.
2
x

Comments on :
(1) If there are several independent variables, each has a truncation error, e.g., O(x2 + t), we say it's first order in time and
second order in space.
(2) The order of a scheme also depends on the local properties of the function, it may be much less than the formal or theoretical
order near sharp gradients. Recall that Taylor series are valid only for smooth and continuous functions.

(3) Typically, we prefer higher-order scheme because is smaller. However, is not necessarily related to accuracy because, for a
given x, a first-order scheme may give more accurate results because the coefficient is smaller. What does tells us is how
the error will change as we change the resolution. Higher order decreases faster when we decreases x.
e.g., = (x)3 versus = x.
Remember that includes , which might be

4u
, etc and it matters.
x 4

(4) Truncation error is cumulative adds up per time step.

(5) Truncation error is usually much larger than machine round-off error. Also for most problems, (space) >> (time), because
solution often evolves smoothly in time but have rapid changes in space. Time step size is often not as large as we wish
because of stability constraint. For meteorological models, we usually want to increase time step size and decrease grid spacing.
General Method for Deriving FD Expressions
Let's write a generalized form of the Taylor series as

mu (x ) m
ui+1 = m
m!
m= 0 x i

Now, suppose we want to use a 3-point stencil at i-1, i, i+1. Then, we can write a generic expression (PDE = FDE + ) as:
u
= aui1 + bui + cui+1 + O( x) m
x

where a, b and c are unknown constants to be determined and m is the order of the approximation. (General rule: If F.D. spans n points,
you can derive an n-1 order F.D. scheme).
Using our Taylor series for ui-1 , ui and ui+1 , we can write

aui 1 + bui + cui+1 =

a (ui x

u ( x ) 2 2u ( x ) 3 3u
+

+ ...)
x
2! x 2
3! x 3

+ bui
u ( x ) 2 2u ( x ) 3 3u
+
+
+ ...)
x
2! x 2
3! x 3
u
( x) 2 2u
( x) 3 3u
= (a + b + c )ui + ( a + c) x + ( a + c )
+
(

a
+
c
)
+ ... (6)
x
2! x 2
3! x 3
c (ui + x

Since we want (6) to have the form of

u
+ O ( x ) m , therefore we set
x

(a+b+c) =0,
(-a+c)x=1
a+c=0.
From them we can find b=0, c = -a = 1/(2x), therefore
u ui+1 ui1
=
+ O ( x) 3
x
2 x

which is a second-order centered difference scheme.

This method can be used in a general manner for symmetric or non-symmetric difference formula. You just pick which points you
want to use.

2.2.2. Polynomial Fitting

This is the second, most general method for generating finite difference expression. Here, we assume that the solution to the PDE can
be approximated by a polynomial, and that the values at the mesh points at exact. We thus differentiate the polynomial to obtain
expression for various derivatives.

10

Assume that u(x) = a x2 + b x + c:

Goal: Find a, b and c. Note that the grid spacing need not be uniform.
Applying the polynomial to those three points gives
ui-1 = a x2 i-1 + b xi-1 + c
ui = a x2 i + b xi + c
ui+1 = a x2 i+1 + b xi+1 + c
Solve for a, b, c, we obtain
( x xi )( x xi+1 )
( x xi1 )( x xi +1 )
u ( x ) = ui1
+ ui

( xi 1 xi )( xi1 xi+1)
(xi xi1 )( xi xi +1 )
( x xi )(x xi1 )
+ ui+1

( xi+1 xi )( xi+1 xi1 )

11

Note: u(xi) = ui for i, i+1, i-1 (verify yourself).

The above formula is often called a Lagrange Interpolation Polynomial and can be generalized to any order.
If the true solution u is a polynomial having the same degree as that used in the interpolation, then the F.D. expression will be
exact as will be the derivatives.
Now, let's compute the derivative

u
x

by differentiating the expression for u(x). First, let us define x xi+1 - xi = xi xi-1 => 2 x =
xi

xi+1 xi-1 , therefore

u ui+1 ui 1
=
x
2 x
2
u ui+1 2ui + ui1
=
x 2
x 2
3u
=0
x 3

because we used a 2nd order polynomial. To obtain FD formulation for third order derivative, we have to use 3rd- or higher- order
polynomial.
The 2nd order case can be interpreted physically as the derivative of derivative, or the difference between two slopes

u 2 u
=
x x x 2
or
ui +1 ui ui ui 1

x
x
x

12

i.e., the difference between the slope centered at i+1/2 and slope centered at i-1/2.

It turns out the a 2nd-order polynomial is usually adequate

Higher-order approximations tend to introduce noises into the solution because the higher-order derivatives may be large,
particularly near sharp gradients.

Beyond 2nd-order, the polynomial method is not identical to the Taylor series method.

2.2.3. Interpolation and Characteristics

Interpolation is at the heart of virtually all numerical techniques. e.g., 1-D advection equation:
u
u
+c
=0
t
x

( c > 0andconstant)

Let's draw a space-time diagram, using discrete points in space and time:

13

Recall that dx/dt = c and du=0 along the characteristic curves.

Goal: Determine uin+1 = f (uin , uin1, uin+1, etc) .
We recognize uin+1 = value of u along characteristic curve at the previous time level, i.e., u*n , i.e.,
uin+1 = u*n .

In general, u* won't coincide with grid points but we know only values at the grid points. So we need interpolation from the grid points
to *.
We are free to use as many points as we want, more points higher order.
Indeed, the interpolation used determines the order of accuracy of the solution.
With simple linear interpolation, we get (see figure)
u =
*

x n x n
ui 1 + 1
ui .
x
x

14

Now, in one time step t, a particle moves a distance x at speed c, i.e.,

x = c t.
Substituting, we have
u*n = uin+1 =

ct n
ct n
ui1 + uin
ui
x
x

or
uin+1 uin
u n uin1
+c i
=0
t
x
This is a familiar first-order upstream method using a forward-in-time scheme!
As we will see later, this scheme is strongly damping, therefore inaccurate. This is because linear interpolation always under-estimate
the peak values in the solution.
This scheme yields exact solution if c = constant and t = x/c. For general problem this condition is hard to meet, however.
If we new include point i+1 in our interpolation, we end up fitting a parabola to the points and obtain
u*n = uin+1 = uin
where =

n
2 n
(ui+1 uin1 ) +
(ui+1 2uin + uin1)
2
2

c t
.
x

This scheme is known as the Lax-Wendroff or the Crowley scheme in 1-D.

Notice that the last term on R.H.S. is an approximation to a second-order derivative which represents a diffusion process. This term is
not present in the original advection equation. It is needed, however, to keep the scheme stable (by damping unstable modes). Without
this term, the scheme becomes forward-in-time and centered-in-space. We will see later in the section of stability analysis that such as

15

scheme is absolutely unstable (i.e., the numerical solution will grow without bound in computer problems, you will quickly run into
floating-point overflow error).

Discussion on Stability Condition / Constraint on Time Step Size

It turns out, as we will show later, that we can only use interpolation not extrapolation at time n for numerical stability.
The criterion of interpolation is that
ct
< 1 , i.e.,
x

|x| < |x|

or particle cannot move more than one grid interval per time step.
This is known as the stability constraint, and we can write it as
t < x/c.
It is a limit on the time step size in terms of x and c. All explicit schemes have stability limitations (all of which are not the same).
The faster signal moves, the smaller is the time step size that can be used. Half x half t. In 3D and for a fixed domain size,
doubling resolution in all three dimensions increases the total computation by a factor (2)4 = 16!

16