<< . .

. 54
( : 59)

. . >>

of c is optimal.

(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
9.3. Finite Volume Methods 383

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.3 Given an arbitrary, but ¬xed, triangle K with diameter hK , prove the
∆p 0,K ¤ |p|1,K
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.

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
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 :=

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
„∈ ’ ,
x = x(„ ) = aij + „ dij νij , ,
and if we introduce the composite function w(„ ) := u(x(„ )) , then we can
µ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
that q is continuous on the interval ’ 2 , 1 , the equation can be solved

dij γij dij 1 1
q(s) exp ’ ds + w ’
w(„ ) = s+
µij µij 2 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 ):

γij dij
w( 1 ) ’ w(’ 1 ) exp
2 2 µij
qij ≈ γij . (9.25)
γij dij
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
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
γij dij
then with the choice rij := R , (9.25) can be written as

uj ’ ui
qij ≈ µij ’ [rij ui + (1 ’ rij ) uj ] γij .
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
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,
the condition (Ah )ii ≥ 0 leads to the requirement
γij dij
’ ¤ 1,
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.
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
¤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
√ 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

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
+ Lu = f in QT ,
u=0 on ST ,
„¦ — {0} ,
u = u0 on

(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
t ∈ (0, T ) ,
X(x, s, t) = c(X(x, s, t), t) ,
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
+ c · ∇u (X(x, s, t), t) .
(x, t) =
‚t ‚t
The particular value
ˆ ‚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 ;
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
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
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.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

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
[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.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

<< . .

. 54
( : 59)

. . >>