(ii) In general, the seminorm |u|k+1 depends on negative powers of µ.

Therefore, if h ’ 0, the convergence in Theorem 9.3 is not uniform

with respect to µ.

Comparing the estimate from Theorem 9.3 for the special case of continuous

linear elements with the estimate (9.13) for the corresponding standard

method given at the end of the introduction, i.e.,

u ’ uh ¤ Ch|u|2 ,

µ

we see that the error of the streamline-di¬usion method is measured in a

stronger norm than the · µ -norm and additionally, that the error bound

is asymptotically better in the interesting case µ < h . A further advan-

tage of the streamline-di¬usion method is to be seen in the fact that its

implementation is not much more di¬cult than that of the standard ¬nite

element method.

However, there are also some disadvantages: Since the error bound in-

volves the H k+1 -seminorm of the solution u, it may depend on negative

powers of µ. Furthermore, there is no general rule to determine the param-

eters δ1 , δ2 . Usually, they are chosen more or less empirically. This may be

a problem when the streamline-di¬usion method is embedded into more

complex programs (for example, for solving nonlinear problems). Finally,

in contrast to the ¬nite volume methods described in the next section, the

property of inverse monotonicity (cf. Theorem 6.19) cannot be proven in

general.

9.3. Finite Volume Methods 383

Exercises

9.1 (a) Given a constant di¬usion coe¬cient µ > 0, rewrite the ordinary

boundary value problem

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

u(0) = u(1) ’ 1 = 0 ,

into an equivalent form but with nonnegative right-hand side and

homogeneous Dirichlet boundary conditions.

(b) Compute the H 2 (0, 1)-seminorm of the solution of the transformed

problem and investigate its dependence on µ.

9.2 Prove the error equation of the streamline-di¬usion method (Corollary

9.2).

9.3 Given an arbitrary, but ¬xed, triangle K with diameter hK , prove the

inequality

cinv

∆p 0,K ¤ |p|1,K

hK

for arbitrary polynomials p ∈ Pk (K), k ∈ N, where the constant cinv > 0

is independent of K and p.

·

9.4 Verify that the streamline-di¬usion norm is indeed a norm.

sd

9.3 Finite Volume Methods

In the convection-dominated situation, the ¬nite volume method intro-

duced in Chapter 6 proves to be a very stable, but not so accurate, method.

One reason for this stability lies in an appropriate asymptotic behaviour of

the weighting function R for large absolute values of its argument.

Namely, if we consider the examples of nonconstant weighting functions

given in Section 6.2.2, we see that

(P4) lim R(z) = 0 , lim R(z) = 1 .

z’’∞ z’∞

In the general case of the model problem (6.5) with k = µ > 0, (P4)

γij dij

’1 the term rij ui + (1 ’ rij )uj in the bilinear

implies that for

µ

γij dij

form bh e¬ectively equals uj , whereas in the case 1 the quantity

µ

ui remains.

In other words, in the case of dominating convection, the approximation

bh evaluates the “information” (uj or ui ) upwind, i.e., just at that node (aj

or ai ) from which “the ¬‚ow is coming”.

384 9. Discretization of Convection-Dominated Problems

This essentially contributes to the stabilization of the method and makes

it possible to prove properties such as global conservativity or inverse mono-

tonicity (cf. Section 6.2.4) without any restrictions on the size of the local

γij dij

P´clet number

e and thus without any restrictions on the ratio of h

µ

and µ. This local P´clet number (note the missing factor 2 in comparison

e

to (9.23)) also takes the direction of the ¬‚ow compared to the edge ai aj

into account.

The Choice of Weighting Parameters

In order to motivate the choice of the weighting parameters in the case of

the Voronoi diagram, we recall the essential step in the derivation of the

¬nite volume method, namely the approximation of the integral

[µij (νij · ∇u) ’ γij u] dσ .

Iij :=

“ij

It ¬rst suggests itself to apply a simple quadrature rule, for example

Iij ≈ qij mij ,

where qij denotes the value of the expression to be integrated at the

point aij of the intersection of the boundary segment “ij with the edge

bounded by the vertices ai and aj (i.e., 2aij = ai + aj ). Next, if this edge

is parametrised according to

11

„∈ ’ ,

x = x(„ ) = aij + „ dij νij , ,

22

and if we introduce the composite function w(„ ) := u(x(„ )) , then we can

write

µij dw

µij (νij · ∇u) ’ γij u = q(0) with q(„ ) := („ ) ’ γij w(„ ) .

dij d„

The relation de¬ning the function q can be interpreted as a linear ordinary

di¬erential equation for the unknown function w : ’ 2 , 1 ’ R . Provided

1

2

that q is continuous on the interval ’ 2 , 1 , the equation can be solved

1

2

exactly:

„

dij γij dij 1 1

q(s) exp ’ ds + w ’

w(„ ) = s+

µij µij 2 2

’1/2

γij dij 1

— exp „+ .

µij 2

Approximating q by a constant qij , we get in the case γij = 0,

qij γij dij 1 1

≈ 1 ’ exp ’ +w ’

w(„ ) „+

γij µij 2 2

γij dij 1

— exp „+ .

µij 2

9.3. Finite Volume Methods 385

In particular,

1 qij γij dij 1 γij dij

≈ 1 ’ exp ’ +w ’

w exp ; (9.24)

2 γij µij 2 µij

that is, the approximation qij of q(0) can be expressed by means of the

values w(± 1 ):

2

γij dij

w( 1 ) ’ w(’ 1 ) exp

2 2 µij

qij ≈ γij . (9.25)

γij dij

’1

exp µij

In the case γij = 0, it immediately follows from the exact solution and the

approximation q ≈ qij that

w( 1 ) ’ w(’ 1 )

qij ≈ µij 2 2

.

dij

Since this is equal to the limit of (9.25) for γij ’ 0, we can exclusively

work with the representation (9.25).

If we de¬ne the weighting function R : R ’ [0, 1] by

1 z

R(z) := 1 ’ 1’ , (9.26)

ez ’ 1

z

γij dij

then with the choice rij := R , (9.25) can be written as

µij

uj ’ ui

qij ≈ µij ’ [rij ui + (1 ’ rij ) uj ] γij .

dij

A simple algebraic manipulation shows that this is exactly the approxima-

tion scheme given in Section 6.2.

The use of the weighting function (9.26) yields a discretization method

that can be interpreted as a generalization of the so-called Il™in“ Allen“

Southwell scheme. However, in order to avoid the comparatively expensive

computation of the function values rij of (9.26), often simpler functions

R : R ’ [0, 1] are used (see Section 6.2.2), which are to some extent

approximations of (9.26) keeping the properties (P1) to (P4).

At the end of this paragraph we will illustrate the importance of the prop-

erties (P1) to (P3), especially for convection-dominated problems. Property

(P2) has been used in the proof of the basic stability estimate (6.20). On

the other hand, we have seen at several places (e.g., in Section 1.4 or in

Chapter 5) that the matrix Ah of the corresponding system of linear al-

gebraic equations should have positive diagonal entries. For example, if

in the di¬erential equation from (9.2) the reaction term disappears, then

properties (P1) and (P3) guarantee that the diagonal entries are at least

nonnegative. This can be seen as follows:

386 9. Discretization of Convection-Dominated Problems

From (6.9) we conclude the following formula:

µij µij γij dij

i ∈ Λ.

(Ah )ii = + γij rij mij = 1+ rij mij ,

dij dij µij

If we replace in property (P3) the number z by ’z, then we get, by property

(P1),

0 ¤ 1 + [1 ’ R(’z)] z = 1 + zR(z) .

Therefore, if the weighting function R satis¬es (P1) and (P3), then we have

that (Ah )ii ≥ 0 for all i ∈ Λ .

The simple choice rij ≡ 1 does not satisfy property (P3). In this case,

2

the condition (Ah )ii ≥ 0 leads to the requirement

γij dij

’ ¤ 1,

2µij

which in the case γij ¤ 0, i.e., for a local ¬‚ow from aj to ai , is a restriction

to the ratio of h and µ, and this is analogous to the condition (9.5) on the

grid P´clet number, where only the sizes of K, c, and h enter.

e

Similarly, it can be shown that property (P3) implies the nonpositivity

of the o¬-diagonal entries of Ah .

An Error Estimate

At the end of this section an error estimate will be cited, which can be

derived similarly to the corresponding estimate of the standard method.

The only special aspect is that the dependence of the occurring quantities

on µ is carefully tracked (see [40]).

Theorem 9.5 Let {Th }h be a regular family of conforming triangula-

tions, all triangles of which are nonobtuse. Furthermore, in addition to the

assumptions on the coe¬cients of the bilinear form (9.9), let f ∈ C 1 („¦).

If the exact solution u of the model problem belongs to H 2 („¦) and if

uh ∈ Vh denotes the approximative solution of the ¬nite volume method

(6.11), where the approximations γij , respectively ri , are chosen according

¯

to (6.7), respectively (6.8), then for su¬ciently small h > 0 the estimate

h

¤C√ [ u h ∈ (0, ¯ ,

u ’ uh + |f |1,∞ ] , h]

µ 2

µ

¯

holds, where both the constant C > 0 and h > 0 do not depend on µ.

In special, but practically not so relevant, cases (for example, if the trian-

gulations are of Friedrichs“Keller type), it is possible to remove the factor

1

√ in the bound above.

µ

Comparing the ¬nite volume method with the streamline-di¬usion

method, we see that the ¬nite volume method is less accurate. However, it

is globally conservative and inverse monotone.

9.4. The Lagrange“Galerkin Method 387

Exercise

9.5 Using an equidistant grid, formulate both the streamline-di¬usion

method and the ¬nite volume method for a one-dimensional model problem

(d = 1, „¦ = (0, 1), r = 0) with constant coe¬cients and compare the re-

sulting discretizations. Based on that comparison, what can be said about

the choice of the parameters in the streamline-di¬usion method?

9.4 The Lagrange“Galerkin Method

In the previous sections, discretization methods for stationary di¬usion-

convection equations were presented. In conjunction with the method of

lines, these methods can also be applied to parabolic problems. However,

since the method of lines decouples spatial and temporal variables, it can-

not be expected that the peculiarities of nonstationary di¬usion-convection

equations are re¬‚ected adequately.

The so-called Lagrange“Galerkin method attempts to bypass this prob-

lem by means of an intermediate change from the Eulerian coordinates

(considered up to now) to the so-called Lagrangian coordinates. The latter

are chosen in such a way that the origin of the coordinate system (i.e., the

position of the observer) is moved with the convective ¬eld, and in the new

coordinates no convection occurs.

To illustrate the basic idea, the following initial-boundary value problem

will be considered, where „¦ ‚ Rd is a bounded domain with Lipschitz

continuous boundary and T > 0:

For given functions f : QT ’ R and u0 : „¦ ’ R , ¬nd a function

u : QT ’ R such that

‚u

+ Lu = f in QT ,

‚t

(9.27)

u=0 on ST ,

„¦ — {0} ,

u = u0 on

where

(Lu) (x, t) := ’∇·(K(x) ∇u(x, t))+c(x, t)·∇u(x, t)+r(x, t)u(x, t), (9.28)

with su¬ciently smooth coe¬cients

K : „¦ ’ Rd,d , c : Q T ’ Rd , r : QT ’ R .

As usual, the di¬erential operators ∇ and ∇· act only with respect to the

spatial variables.

The new coordinate system is obtained by solving the following

parameter-dependent auxiliary problem:

388 9. Discretization of Convection-Dominated Problems

Given (x, s) ∈ QT , ¬nd a vector ¬eld X : „¦ — [0, T ]2 ’ Rd such that

d

t ∈ (0, T ) ,

X(x, s, t) = c(X(x, s, t), t) ,

(9.29)

dt

X(x, s, s) = x .

The trajectories X(x, s, ·) are called characteristics (through (x, s)). If c

is continuous on QT and, for ¬xed t ∈ [0, T ], Lipschitz continuous with

respect to the ¬rst argument on „¦, then there exists a unique solution

X = X(x, s, t). Denoting by u the su¬ciently smooth solution of (9.27)

and setting

u(x, t) := u(X(x, s, t), t) for ¬xed s ∈ [0, T ] ,

ˆ

then the chain rule implies that

‚u

ˆ ‚u

+ c · ∇u (X(x, s, t), t) .

(x, t) =

‚t ‚t

The particular value

‚u

ˆ ‚u

(x, s) + c(x, s) · ∇u(x, s)

(x, s) =

‚t ‚t

is called the material derivative of u at (x, s). Thus the di¬erential equation

reads as

‚u ˆ

’ ∇ · (K ∇u) + ru = f ;

‚t

i.e., it is formally free of any convective terms.

Now the equation will be semidiscretized by means of the horizontal

method of lines. A typical way is to approximate the time derivative by

backward di¬erence quotients. So let an equidistant partition of the time

interval (0, T ) with step size „ := T /N, N ∈ N (provided that T < ∞), be

given.

Tracking the characteristics backwards in time, in the strip „¦ —

[tn , tn+1 ), n ∈ {0, 1, . . . , N ’ 1}, with x = X(x, tn+1 , tn+1 ) the following

approximation results:

‚u ˆ 1 1

≈ [ˆ(x, tn+1 ) ’ u(x, tn )] = [u(x, tn+1 ) ’ u(X(x, tn+1 , tn ), tn )] .

u ˆ

‚t „ „

Further, if Vh denotes a ¬nite-dimensional subspace of V in which we want

to ¬nd the approximations to u(·, tn ), the method reads as follows:

Given u0h ∈ Vh , ¬nd an element U n+1 ∈ Vh , n ∈ {0, . . . , N ’ 1}, such

that

1

U n+1 ’ U n (X(·, tn+1 , tn )), vh 0

„

+ K∇U n+1 · ∇vh , 1 0 + r(·, tn+1 )U n+1 , vh 0 = f (·, tn+1 ), vh 0

for all vh ∈ Vh ,

U 0 = u0h .

(9.30)

9.4. Lagrange“Galerkin Method 389

A possible extension of the method is to use time-dependent subspaces;

that is, given a sequence of subspaces Vh ‚ V, n ∈ {0, . . . , N }, the

n

approximations U n to u(·, tn ) are chosen from Vh .n

So the basic idea of the Lagrange“Galerkin method, namely, the elimi-

nation of convective terms by means of an appropriate transformation of

coordinates, allows the application of standard discretization methods and

makes the method attractive for situations where convection is dominating.

In fact, there exists a whole variety of papers dealing with error esti-

mates for the method in the convection-dominated case, but often under

the condition that the system (9.29) is integrated exactly.

In practice, the exact integration is impossible, and the system (9.29) has

to be solved numerically (cf. [61]). This may lead to stability problems, so

there is still a considerable need in the theoretical foundation of Lagrange“

Galerkin methods.

Only recently it has been possible for a model situation to prove order

of convergence estimates uniformly in the P´clet number for (9.30) (see

e

[43]). The key is the consequent use of Lagrangian coordinates, reveal-

ing that (9.30) is just the application of the implicit Euler method to an

equation arising from a transformation by characteristics de¬ned piecewise

backward in time. This equation is a pure di¬usion problem, but with a

coe¬cient re¬‚ecting the transformation. In conjunction with the backward

Euler method this is not visible in the elliptic part to be discretized. Thus

tracking the characteristics backward in time turns out to be important.

A

Appendices

A.1 Notation

C set of complex numbers

N set of natural numbers

N0 := N ∪ {0}

Q set of rational numbers

R set of real numbers

R+ set of positive real numbers

Z set of integers

z real part of the complex number z

z imaginary part of the complex number z

transpose of the vector x ∈ Rd , d ∈ N

T

x

1/p

d