resulting discrete problems are too expensive or even untreatable.

So one can ask whether the theory is insu¬cient or not. The following

example will show that this is not necessarily the case.

Example 9.1 Given a constant di¬usion coe¬cient k > 0, consider the

boundary value problem

(’ku + u) = 0 in „¦ := (0, 1) ,

u(0) = u(1) ’ 1 = 0 .

Its solution is

1 ’ exp (x/k)

u(x) = .

1 ’ exp (1/k)

A rough sketch of the graph (Figure 9.1) shows that this function has a

signi¬cant boundary layer at the right boundary of the interval even for

the comparatively small global P´clet number Pe = 100. In the larger

e

subinterval (about (0, 0.95)) it is very smooth (nearly constant), whereas

in the remaining small subinterval (about (0.95, 1)) the absolute value of

its ¬rst derivative is large.

Given an equidistant grid of width h = 1/(M + 1), M ∈ N, a dis-

cretization by means of symmetric di¬erence quotients yields the di¬erence

370 9. Discretization of Convection-Dominated Problems

1

0.8

0.6

0.4

0.2

0

0 0.2 0.4 0.6 0.8 1

Figure 9.1. Solution for k = 0.01.

equations

ui’1 ’ 2ui + ui+1 ui+1 ’ ui’1

’k i ∈ {1, . . . , M } =: Λ ,

+ = 0,

h2 2h

u0 = uM+1 ’ 1 = 0 .

Collecting the coe¬cients and multiplying the result by 2h, we arrive at

2k 4k 2k

’ ’ 1 ui’1 + ui + ’ i ∈ Λ.

+ 1 ui+1 = 0 ,

h h h

If we make the ansatz ui = »i , the di¬erence equations can be solved

exactly:

i

1’ 2k+h

2k’h

ui = .

M+1

1’ 2k+h

2k’h

In the case 2k < h, which is by no means unrealistic (e.g., for the typical

value k = 10’7 ), the numerical solution considerably oscillates, in contrast

to the behaviour of the exact solution u. These oscillations do not disappear

until h < 2k is reached, but this condition is very restrictive for small values

of k.

But even if the condition h < 2k is satis¬ed, undesirable e¬ects can

be observed. For example, in the special case h = k we have at the node

aM = M h that

1 ’ exp (M h/k) 1 ’ exp (M ) exp (’M ) ’ 1

u(aM ) = = =

1 ’ exp (1/k) 1 ’ exp (M + 1) exp (’M ) ’ exp (1)

’ exp (’1) for h ’ 0 ,

9.1. Standard Methods 371

whereas the numerical solution at this point asymptotically behaves like

(note that » = (2k + h)/(2k ’ h) = 3)

»’M ’ 1

1 ’ »M 1 1

’= for h ’ 0 .

uM = = ’M

1’» ’»

M+1 » » 3

So the numerical solution does not converge to the exact solution at the

node aM .

Again this is no contradiction to possible convergence results for the ¬nite

di¬erence method in the discrete maximum norm, since now the di¬usion

coe¬cient is not ¬xed, but rather the discretization is to be viewed as

belonging to the limit case k = 0, with an arti¬cial di¬usion part in the

discretization (see (9.8) below).

Finite Di¬erence Methods with Symmetric and One-Sided Dif-

ference Quotients

The oscillations in Example 9.1 show that in this case no comparison prin-

ciple as in Corollary 1.13 is valid. Such a comparison principle, or more

strongly a maximum principle, will lead to nonnegative solutions in the

case of nonnegative right-hand side and Dirichlet data. This avoids for a

homogeneous right-hand side an undershooting, as observed in Example 9.1,

i.e., negative solution values in this case, and also an overshooting, i.e., so-

lution values larger than the maximum of the Dirichlet data, provided that

condition (1.32) (6)* holds.

In the following we will examine how the convective part in¬‚uences the

matrix properties (1.32) and thus the validity of a maximum or comparison

principle and also conclude a ¬rst simple remedy.

We consider the model problem (9.2), for simplicity on a rectangle

„¦ = (0, a) — (0, b), with constant, scalar K = kI and equipped with an

equidistant grid „¦h . To maintain the order of consistency 2 of a spatial

discretization of ’∇ · (K∇u) = ’k∆u by the ¬ve-point stencil, the use of

the symmetric di¬erence quotient for the discretization of

(c · ∇u)(x) = c1 (x)‚1 u(x) + c2 (x)‚2 u(x)

suggests itself, i.e., for a grid point x ∈ „¦h ,

1

c1 (x)‚1 u(x) ∼ c1 (x) (ui+1,j ’ ui’1,j ), (9.3)

2h

and similarly for c2 (x)‚2 u(x) (cf. (1.7) for the notation). This leads to

˜

the following entries of the system matrix Ah , for example in a rowwise

372 9. Discretization of Convection-Dominated Problems

numbering (compare (1.13)):

c1 (x) k

’ ’ 2;

left secondary diagonal:

2h h

c1 (x) k

’ 2;

right secondary diagonal: +

2h h

c2 (x) k

’ ’ 2;

l + 1 positions to the left:

2h h

c2 (x) k

’ 2;

l + 1 positions to the right: +

2h h

4k

diagonal: .

h2

Condition (1.32) (1) and (1.32) (6)* obviously hold.

We check the conditions su¬cient for a comparison principle (Corol-

lary 1.13). To satisfy condition (1.32) (2) we require

|c1 (x)|

k

’ + < 0,

h2 2h

|c2 (x)|

k

’ 2+ < 0.

h 2h

Denoting the grid P´clet number by

e

c ∞h

Peh := , (9.4)

2k

the above conditions are satis¬ed if

Peh < 1 (9.5)

is satis¬ed. Under this assumption also the conditions (1.32) (5) and (7)

are satis¬ed, and thus also (3), i.e., (9.5), is su¬cient for the validity of a

comparison principle. In Example 9.1 this is just the condition h < 2k.

The grid P´clet number is obviously related to the global P´clet number

e e

from (9.1) by

h

Peh = Pe .

2 diam(„¦)

The requirement (9.4) can always be met by choosing h su¬ciently small,

but for large Pe this may be a severe requirement, necessary for the sake of

stability of the method, whereas for the accuracy desired a larger step size

may be su¬cient. A simple remedy to ensure condition (1.32) (2) is to use

a one-sided (upwind) discretization of c1 ‚1 u and c2 ‚2 u, which is selected

against the stream direction de¬ned by c1 and c2 , respectively:

1

For c1 (x) ≥ 0 : c1 (x)‚1 u(x) ∼ c1 (x) (ui,j ’ ui’1,j ) , (9.6)

h

1

c1 (x)‚1 u(x) ∼ c1 (x) (ui+1,j ’ ui,j ) ,

for c1 (x) < 0 :

h

9.1. Standard Methods 373

and analogously for c2 ‚2 u.

Due to this choice there are only additional nonnegative addends to the

diagonal position and nonpositive ones to the o¬-diagonal positions com-

pared to the ¬ve-point stencil or another discretization of a di¬usive part.

Thus all properties (1.32) (1)“(7), (4)*, and (6)* remains una¬ected; i.e.,

the upwind discretization satis¬es all qualitative properties of Section 1.4

from the inverse monotonicity to the strong maximum principle, without

any restrictions to the local P´clet number.

e

The drawback lies in the reduced accuracy, since the one-sided di¬erence

quotient has only order of consistency 1. In Section 9.3 we will develop

more re¬ned upwind discretizations.

Due to

u(x, y) ’ u(x ’ h, y)

c1 = (9.7)

h

c1 h ’u(x ’ h, y) + 2u(x, y) ’ u(x + h, y) u(x + h, y) ’ u(x ’ h, y)

+ c1 ,

h2

2 2h

and analogously for the forward di¬erence quotient, the upwind dis-

cretization can be perceived as a discretization with symmetric di¬erence

quotients if a step-size-dependent di¬usive part, also discretized with ‚ ’ ‚ + ,

is added with the di¬usion coe¬cient

|c1 (x)|

h 0

Kh (x) := . (9.8)

|c2 (x)|

0

2

Therefore, one also speaks of adding arti¬cial di¬usion (or viscosity). The

disadvantage of this full upwind method is that it recognizes the ¬‚ow di-

rection only if the ¬‚ow is aligned to one of the coordinate axes. This will

be improved in Section 9.2.

Error Estimates for the Standard Finite Element Method

In order to demonstrate the theoretical de¬ciencies, we will again reproduce

the way for obtaining standard error estimates for a model problem. So let

K(x) ≡ µI with a constant coe¬cient µ > 0, c ∈ C 1 („¦, Rd ), r ∈ C(„¦),

f ∈ L2 („¦). Furthermore, assume that the following inequality is valid in

„¦, where r0 > 0 is a constant: r ’ 1 ∇ · c ≥ r0 .

2

Then the bilinear form a : V — V ’ R, V := H0 („¦), corresponding to

1

the boundary value problem (9.2), reads as (cf. (3.23))

[µ∇u · ∇v + c · ∇u v + r uv ] dx , u, v ∈ V .

a(u, v) := (9.9)

„¦

To get an ellipticity estimate of a, we set u = v ∈ V in (9.9) and take the

relation 2v(c · ∇v) = c · ∇v 2 into account. Then, by partial integration of

the middle term, we obtain

µ|v|2 + c · ∇v, v

a(v, v) = + rv, v

1 0 0

374 9. Discretization of Convection-Dominated Problems

1 1

µ|v|2 ’ ∇ · c, v 2 = µ|v|2 + r ’ ∇ · c, v 2

= + rv, v .

1 1

0

2 2

0 0

Introducing the so-called µ-weighted H 1 -norm by

2 1/2

:= µ|v|2 + v

v , (9.10)

µ 1 0

we immediately arrive at the estimate

a(v, v) ≥ µ|v|2 + r0 v ≥ ± v 2,

2

˜ (9.11)

1 0 µ

where ± := min{1, r0 } does not depend on µ.

˜

Due to c · ∇u = ∇ · (cu) ’ (∇ · c)u, partial integration yields for arbitrary

u, v ∈ V the identity

c · ∇u, v = ’ u, c · ∇v ’ (∇ · c)u, v .

0 0 0

So we get the continuity estimate

|a(u, v)| ¤ µ|u|1 |v|1 + c 0,∞ u 0|v|1 + (|c|1,∞ + r 0,∞ ) u 0 v 0

√ √

¤ ( µ |u|1 + u 0 ) {( µ + c 0,∞ ) |v|1 + (|c|1,∞ + r 0,∞ ) v 0 }

¤ M u µ v 1,

˜

(9.12)

˜ := 2 max{√µ + c 0,∞ , |c|1,∞ + r 0,∞ }.

where M

Since we are interested in the case of small di¬usion µ > 0 and present

˜

convection (i.e., c 0,∞ > 0), the continuity constant M can be bounded

independent of µ. It is not very surprising that the obtained continuity

estimate is nonsymmetric, since also the di¬erential expression L behaves

like that. Passing over to a symmetric estimate results in the following

relation:

˜

M

√ u µ v µ.

|a(u, v)| ¤

µ

Now, if Vh ‚ V denotes a ¬nite element space, we can argue as in the proof

of C´a™s lemma (Theorem 2.17) and get an error estimate for the corre-

e

sponding ¬nite element solution uh ∈ Vh . To do this, the nonsymmetric

continuity estimate (9.12) is su¬cient. Indeed, for arbitrary vh ∈ Vh , we

have

± u ’ uh ¤ a(u ’ uh, u ’ uh) = a(u ’ uh , u ’ vh ) ¤ M u ’ uh u ’ vh

˜

2

˜ .

µ 1

µ

Thus

˜

M

u ’ uh ¤ inf u ’ vh 1.

µ

± vh ∈Vh

˜

˜˜

Here the constant M /± does not depend on µ, h, and u. This estimate

is weaker than the standard estimate, because the µ-weighted H 1 -norm is

weaker than the H 1 -norm. Moreover, the error of the best approximation is

9.2. The Streamline-Di¬usion Method 375

not independent of µ, in general. For example, if we apply continuous, piece-

wise linear elements, then, under the additional assumption u ∈ H 2 („¦),

Theorem 3.29 yields the estimate

u ’ vh ¤ u ’ Ih (u) ¤ Ch|u|2 ,

inf 1 1

vh ∈Vh

where the constant C > 0 does not depend on µ, h, and u. So, we ¬nally

arrive at the relation

u ’ uh ¤ Ch|u|2 . (9.13)

µ

However, the H 2 -seminorm of the solution u depends on µ in a disadvan-

tageous manner; for example, it may be (cf. also [27, Lemma III.1.18])

that

|u|2 = O(µ’3/2 ) (µ ’ 0) .

This result is sharp, since for examples of boundary value problems for

ordinary linear di¬erential equations the error of the best approximation

already exhibits this asymptotic behaviour.

So the practical as well as the theoretical problems mentioned above

indicate the necessity to use special numerical methods for solving

convection-dominated equations. In the next sections, a small collection

of these methods will be depicted.

9.2 The Streamline-Di¬usion Method

The streamline-di¬usion method is the prevalent method in the numerical

treatment of stationary convection-dominated problems. The basic idea

is due to Brooks and Hughes [49], who called the method the streamline

upwind Petrov“Galerkin method (SUPG method).

We describe the idea of the method for a special case of boundary value

problem (9.2) under consideration. Let the domain „¦ ‚ Rd be a bounded

polyhedron. We consider the same model as in the preceding section, that

is, K(x) ≡ µI with a constant coe¬cient µ > 0, c ∈ C 1 („¦, Rd ), r ∈ C(„¦),

f ∈ L2 („¦). We also assume that the inequality r ’ 1 ∇ · c ≥ r0 is valid in „¦,

2

where r0 > 0 is a constant. Then the variational formulation of (9.2) reads

as follows:

Find u ∈ V such that

for all v ∈ V ,

a(u, v) = f, v (9.14)

0

where a is the bilinear form (9.9).

Given a regular family of triangulations {Th }, let Vh ‚ V denote the set

of continuous functions that are piecewise polynomial of degree k ∈ N and

satisfy the boundary conditions, i.e.,

Vh := vh ∈ V vh |K ∈ Pk (K) for all K ∈ Th . (9.15)

376 9. Discretization of Convection-Dominated Problems

If in addition the solution u ∈ V of (9.14) belongs to the space H k+1 („¦),

we have, by (3.87), the following error estimate for the interpolant Ih (u):

u ’ Ih (u) ¤ cint hk+1’l |u|k+1,K (9.16)

l,K K

for 0 ¤ l ¤ k+1 and all K ∈ Th . Since the spaces Vh are of ¬nite dimension,

a so-called inverse inequality can be proven (cf. Theorem 3.43, (2) and

Exercise 9.3):

cinv

∆vh 0,K ¤ |vh |1,K (9.17)

hK

for all vh ∈ Vh and all K ∈ Th . Here it is important that the constants

cint , cinv > 0 from (9.16) and (9.17), respectively, do not depend on u or vh

and on the particular elements K ∈ Th .

The basic idea of the streamline-di¬usion method consists in the addition

of suitably weighted residuals to the variational formulation (9.14). Because

of the assumption u ∈ H k+1 („¦), k ∈ N, the di¬erential equation can be

interpreted as an equation in L2 („¦). In particular, it is valid on any element

K ∈ Th in the sense of L2 (K), i.e.,

’µ∆u + c · ∇u + ru = f almost everywhere in K and for all K ∈ Th .

Next we take an elementwise de¬ned mapping „ : Vh ’ L2 („¦) and multiply

the local di¬erential equation in L2 (K) by the restriction of „ (vh ) to K.

Scaling by a parameter δK ∈ R and summing the results over all elements

K ∈ Th , we obtain

δK ’µ∆u + c · ∇u + ru, „ (vh ) = δK f, „ (vh ) .

0,K 0,K

K∈Th K∈Th