<< . .

. 52
( : 59)



. . >>

may happen that the theoretically required step sizes are so small that the
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

<< . .

. 52
( : 59)



. . >>