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

If we add this relation to equation (9.14) restricted to Vh , we see that the

weak solution u ∈ V © H k+1 („¦) satis¬es the following variational equation:

for all vh ∈ Vh ,

ah (u, vh ) = f, vh h

where

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

ah (u, vh ) := a(u, vh ) + ,

0,K

K∈Th

f, vh := f, v + δK f, „ (vh ) .

h 0 0,K

K∈Th

Then the corresponding discretization reads as follows:

Find uh ∈ Vh such that

for all vh ∈ Vh .

ah (uh , vh ) = f, vh (9.18)

h

Corollary 9.2 Suppose the problems (9.14) and (9.18) have a solution u ∈

V © H k+1 („¦) and uh ∈ Vh , respectively. Then the following error equation

is valid:

ah (u ’ uh , vh ) = 0 for all vh ∈ Vh . (9.19)

9.2. Streamline-Di¬usion Method 377

In the streamline-di¬usion method (sdFEM), the mapping „ used in (9.18)

is chosen as „ (vh ) := c · ∇vh .

Without going into details, we mention that a further option is to set

„ (vh ) := ’µ∆vh + c · ∇vh + rvh . This results in the so-called Galerkin/least

squares“FEM (GLSFEM) [54].

Especially with regard to the extension of the method to other ¬nite ele-

ment spaces, the discussion of how to choose „ and δK is not yet complete.

Interpretation of the Additional Term in the Case of Linear Ele-

ments

If the ¬nite element spaces Vh are formed by piecewise linear functions (i.e.,

in the above de¬nition (9.15) of Vh we have k = 1), we get ∆vh |K = 0 for

all K ∈ Th . If in addition there is no reactive term (i.e., r = 0), the discrete

bilinear form is

µ∇uh ·∇vh dx+ c · ∇uh , vh 0 + δK c · ∇uh , c · ∇vh

ah (uh , vh ) = .

0,K

„¦ K∈Th

Since the scalar product appearing in the sum can be rewritten as

c · ∇uh , c · ∇vh 0,K = K (ccT ∇uh ) · ∇vh dx , we obtain the following

equivalent representation:

(µI + δK ccT )∇uh · ∇vh dx + c · ∇uh , vh

ah (uh , vh ) = .

0

K

K∈Th

This shows that the additional term introduces an element-dependent extra

di¬usion in the direction of the convective ¬eld c (cf. also Exercise 0.3),

which motivates the name of the method. In this respect, the streamline-

di¬usion method can be understood as an improved version of the full

upwind method, as seen, for example, in (9.6).

Analysis of the Streamline-Di¬usion Method

To start the analysis of stability and convergence properties of the

streamline-di¬usion method, we consider the term ah (vh , vh ) for arbitrary

vh ∈ Vh .

As in Section 3.2.1, the structure of the discrete bilinear form ah allows

us to derive the estimate

ah (vh , vh ) ≥ µ|vh |2 +r0 vh ’µ∆vh + c · ∇vh + rvh , c · ∇vh

2

0+ δK .

1 0,K

K∈Th

Furthermore, neglecting for a moment the second term in the sum and

using the elementary inequality ab ¤ a2 + b2 /4 for arbitrary a, b ∈ R, we

378 9. Discretization of Convection-Dominated Problems

get

δK ’µ∆vh + rvh , c · ∇vh 0,K

K∈Th

¤ ’µ |δK | ∆vh , |δK | c · ∇vh

0,K

K∈Th

|δK | rvh , |δK | c · ∇vh

+

0,K

¤ µ2 |δK | ∆vh + |δK | r

2 2 2

vh

0,K 0,∞,K 0,K

K∈Th

|δK |

c · ∇vh 2

+ .

0,K

2

By means of the inverse inequality (9.17) it follows that

c2

δK ’µ∆vh + rvh , c · ∇vh ¤ µ2 |δK | 2 |vh |1,K

inv 2

0,K

hK

K∈Th K∈Th

|δK |

+|δK | r c · ∇vh 2

2 2

vh + .

0,K

0,∞,K 0,K

2

Putting things together, we obtain

c2

ah (vh , vh ) ≥ µ ’ µ |δK | inv |vh |2

2 2

1,K vh 0,K

h2

K

K∈Th

|δK |

+ r0 ’ |δK | r + δK ’ c · ∇vh

2 2

.

0,∞,K 0,K

2

The choice

h2

1 r0

¤ min K

0 < δK 2 , r2 (9.20)

2 µcinv 0,∞,K

leads to

µ r0 1

ah (vh , vh ) ≥ |vh |2 + δK c · ∇vh

2 2

vh + .

1 0 0,K

2 2 2

K∈Th

Therefore, if the so-called streamline-di¬usion norm is de¬ned by

1/2

δK c · ∇v v∈V ,

µ|v|2 + r0 v 2 2

v := + ,

sd 1 0 0,K

K∈Th

then the choice (9.20) implies the estimate

1

¤ ah (vh , vh ) for all vh ∈ Vh .

2

vh (9.21)

sd

2

9.2. Streamline-Di¬usion Method 379

Obviously, the streamline-di¬usion norm · sd is stronger than the

µ-weighted H 1 -norm (9.10); i.e.,

√

min{1, r0 } v µ ¤ v sd for all v ∈ V .

Now an error estimate can be proven. Since estimate (9.21) holds only on

the ¬nite element spaces Vh , we consider ¬rst the norm of Ih (u) ’ uh ∈ Vh

and make use of the error equation (9.19):

1

Ih (u) ’ uh ¤ ah (Ih (u) ’ uh , Ih (u) ’ uh ) = ah (Ih (u) ’ u, Ih (u) ’ uh ) .

2

sd

2

In particular, under the assumption u ∈ V © H k+1 („¦) the following three

estimates are valid:

√

∇(Ih (u) ’ u) · ∇(Ih (u) ’ uh ) dx ¤ µ |Ih (u) ’ u|1 Ih (u) ’ uh

µ sd

„¦

√

¤ cint µ hk |u|k+1 Ih (u) ’ uh ,

sd

[c · ∇(Ih (u) ’ u) + r(Ih (u) ’ u)](Ih (u) ’ uh ) dx

„¦

(r ’ ∇ · c)(Ih (u) ’ u)(Ih (u) ’ uh ) dx

=

„¦

’ (Ih (u) ’ u) c · ∇(Ih (u) ’ uh ) dx

„¦

¤ r’∇·c Ih (u) ’ u Ih (u) ’ uh

0,∞ 0 0

+ Ih (u) ’ u c · ∇(Ih (u) ’ uh )

0 0

±

1/2

¤ Ih (u) ’ u 2

C

0,K

K∈Th

1/2

’1

δK Ih (u) ’ u Ih (u) ’ uh

2

+ sd

0,K

K∈Th

1/2

’1

¤ 1 + δK h2 |u|2 Ih (u) ’ uh

Chk ,

sd

K k+1,K

K∈Th

and

δK ’µ∆(Ih (u) ’ u) + c · ∇(Ih (u) ’ u)

K∈Th

+ r(Ih (u) ’ u), c · ∇(Ih (u) ’ uh ) (9.22)

0,K

380 9. Discretization of Convection-Dominated Problems

¤ δK µhk’1 + c k+1

k

cint 0,∞,K hK +r 0,∞,K hK

K

K∈Th

— |u|k+1,K δK c · ∇(Ih (u) ’ uh ) 0,K

1/2

2

¤C |u|2 Ih (u) ’ uh

δK µhk’1 + hk + hk+1 sd .

K k+1,K

K K

K∈Th

Condition (9.20), which was already required for estimate (9.21), implies

that

h2

¤ 2K ,

µδK

cinv

and so the application to the ¬rst term of the last bound leads to

δK ’ µ∆(Ih (u) ’ u) + c · ∇(Ih (u) ’ u)

K∈Th

+ r(Ih (u) ’ u), c · ∇(Ih (u) ’ uh ) 0,K

1/2

¤ Ch δK ] |u|2 Ih (u) ’ uh

k

[µ + .

sd

k+1,K

K∈Th

Collecting the estimates and dividing by Ih (u) ’ uh sd , we obtain the

relation

1/2

h2

Ih (u) ’ uh ¤ Chk µ + K + h2 + δK |u|2 .

sd K k+1,K

δK

K∈Th

Finally, the terms in the square brackets will be equilibrated with the help

of condition (9.20). We rewrite the µ-dependent term in this condition as

h2 2

K

=2 PeK hK

2

µcinv cinv c ∞,K

with

c ∞,K hK

PeK := . (9.23)

2µ

This local P´clet number is a re¬nement of the de¬nition (9.4).

e

The following distinctions concerning PeK are convenient:

PeK ¤ 1 and PeK > 1 .

In the ¬rst case, we choose

h2 2

= δ1 K ,

δK = δ0 PeK hK δ0 = δ1 ,

µ c ∞,K

9.2. Streamline-Di¬usion Method 381

with appropriate constants δ0 > 0 and δ1 > 0, respectively, which are

independent of K and µ. Then we have

h2 1 2PeK

hK ¤ C(µ + hK ) ,

+ h2 + δK = µ + h2 + δ1

K

µ+ 1+

K K

δK δ1 c 0,∞,K

where C > 0 is independent of K and µ. In the second case, it is su¬cient to

choose δK = δ2 hK with an appropriate constant δ2 > 0 that is independent

of K and µ. Then

δ2 c 0,∞,K h2

δ2 K

δK = PeK hK =

PeK 2PeK µ

and

h2 1

+ δ2 hK + h2 ¤ C(µ + hK ) ,

µ + K + h2 + δK = µ +

K K

δK δ2

with C > 0 independent of K and µ. Note that in both cases the constants

can be chosen su¬ciently small, independent of PeK , that the condi-

tion (9.20) is satis¬ed. Now we are prepared to prove the following error

estimate.

Theorem 9.3 Let the parameters δK be given by

±

h2

δ1 K , PeK ¤ 1 ,

δK =

δ hµ , Pe > 1 ,

2K K

where δ1 , δ2 > 0 do not depend on K and µ and are chosen such that

condition (9.20) is satis¬ed. If the weak solution u of (9.14) belongs to

H k+1 („¦), then

√

√

u ’ uh sd ¤ C µ + h hk |u|k+1 ,

where the constant C > 0 is independent of µ, h, and u.

Proof: By the triangle inequality, we get

u ’ uh ¤ u ’ Ih (u) + Ih (u) ’ uh sd .

sd sd

An estimate of the second addend is already known. To deal with the ¬rst

term, the estimates of the interpolation error (9.16) are used directly:

u ’ Ih (u) 2

sd

µ|u ’ Ih (u)|2 + r0 u ’ Ih (u) δK c · ∇(u ’ Ih (u))

2 2

= +

1 0 0,K

K∈Th

2(k+1)

¤ |u|2

c2 µh2k + r0 hK 2 2k

+ δK c 0,∞,K hK

int K k+1,K

K∈Th

382 9. Discretization of Convection-Dominated Problems

¤ µ + h2 + δK |u|2

k+1,K ¤ C(µ + h)hK |u|k+1 .

Ch2k 2k 2

K K

K∈Th

2

Remark 9.4 (i) In the case of large local P´clet numbers, we have

e

µ¤ 1 c ∞,K hK and thus

2

1/2

u ’ uh hK c · ∇(u ’ uh ) ¤ Chk+1/2 |u|k+1 .

2

+ δ2

0 0,K

K∈Th

So the L2 -error of the solution is not optimal in comparison with the

estimate of the interpolation error

u ’ Ih (u) ¤ Chk+1 |u|k+1 ,

0

whereas the L2 -error of the directional derivative of u in the direction

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

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

|x|p j=1 |xj | , x = (x1 , . . . , xd )T ∈ Rd , d ∈ N, p ∈ [1, ∞)

p

:=

|x|∞ := maxj=1,...,d |xj | maximum norm of the vector x ∈ Rd , d ∈ N

|x| := |x|2 Euclidean norm of the vector x ∈ Rd , d ∈ N

d

x·y := xT y = j=1 xj yj scalar product of the vectors x, y ∈ Rd

:= y T Ax = y · Ax energy product of the vectors x, y ∈ Rd w.r.t.

x, y A

a symmetric, positive de¬nite matrix A

|±| := |±|1 order (or length) of the multi-index ± ∈ Nd , d ∈ N 0

I identity matrix or identity operator

jth unit vector in Rm , j = 1, . . . , m

ej

diagonal matrix in Rm,m with diagonal

diag(»i ) = diag(»1 , . . . , »m )

entries »1 , . . . , »m ∈ C

A.1. Notation 391

AT transpose of the matrix A

A’T transpose of the inverse matrix A’1

det A determinant of the square matrix A

»min (A) minimum eigenvalue of a matrix A with real eigenvalues

»max (A) maximum eigenvalue of a matrix A with real eigenvalues

σ(A) set of eigenvalues (spectrum) of the square matrix A

(A) spectral radius of the square matrix A

m(A) bandwidth of the symmetric matrix A

Env (A) hull of the square matrix A

p(A) pro¬le of the square matrix A

:= {x : x ’ x0 < } open ball in a normed space

B (x0 )

:= {x : x ’ x0 ¤ } closed ball in a normed space

B (x0 )

diameter of the set G ‚ Rd

diam (G)

|G|n n-dimensional (Lebesgue) measure of the G ‚ Rn , n ∈ {1, . . . , d}

|G| := |G|d d-dimensional (Lebesgue) measure of the set G ‚ Rd

vol (G) length (d = 1), area (d = 2), volume (d = 3) of “geometric bodies”

G ‚ Rd

int G interior of the set G

‚G boundary of the set G

G closure of the set G

span G linear hull of the set G

conv G convex hull of the set G

|G| cardinal number of the discrete set G

outer unit normal w.r.t. the set G ‚ Rd

ν

domain of Rd , d ∈ N

„¦

:= ‚„¦ boundary of the domain „¦ ‚ Rd

“

supp • support of the function •

’1

f inverse of the mapping f

f [G] image of the set G under the mapping f

’1

f [G] preimage of the set G under the mapping f

f |K restriction of f : G ’ R to a subset K ‚ G

vX norm of the element v of the normed space X

dim X dimension of the ¬nite-dimensional linear space X

L[X, Y ] set of linear, continuous operators acting from the normed space

X in the normed space Y

:= L[X, R] dual space of the real normed space X

X

O(·), o(·) Landau symbols of asymptotic analysis

(i, j ∈ N0 ) Kronecker symbol, i.e., δii = 1 and δij = 0 if i = j

δij

Di¬erential expressions

(l ∈ N) symbol for the partial derivative w.r.t. the lth variable

‚l

(t ∈ R) symbol for the partial derivative w.r.t. the variable t

‚t

(± ∈ Nd multi-index) ±th partial derivative

‚± 0

∇ := (‚1 , . . . , ‚d )T Nabla operator (symbolic vector)

392 A. Appendices

∆ Laplace operator

:= µ · ∇ directional derivative w.r.t. the vector µ

‚µ

:= ‚¦ := (‚j ¦i )m

D¦ Jacobi matrix or functional matrix

i,j=1

‚x

of a di¬erentiable mapping ¦ : Rm ’ Rm

Coe¬cients in di¬erential expressions

K di¬usion coe¬cient (a square matrix function)

c convection coe¬cient (a vector function)

r reaction coe¬cient

Discretization methods

Vh ansatz space

Xh extended ansatz space without any homogeneous Dirichlet

boundary conditions

ah approximated bilinear form

bh approximated linear form

Function spaces (see also Appendix A.5)

Pk (G) set of polynomials of maximum degree k on G ‚ Rd

C(G) = C 0 (G) set of continuous functions on G

(l ∈ N) set of l-times continuously di¬erentiable functions on G

C l (G)

∞

C (G) set of in¬nitely often continuously di¬erentiable functions on G

0

C(G) = C (G) set of bounded and uniformly continuous functions on G

(l ∈ N) set of functions with bounded and uniformly continuous

C l (G)

derivatives up to the order l on G

∞

C (G) set of functions, all partial derivatives of which are bounded and

uniformly continuous on G

0

C0 (G) = C0 (G) set of continuous functions on G with compact support

(l ∈ N) set of l-times continuously di¬erentiable functions on G

l

C0 (G)

with compact support

∞

C0 (G) set of in¬nitely often continuously di¬erentiable functions on G

with compact support

(p ∈ [1, ∞)) set of Lebesgue-measurable functions whose pth power

p

L (G)

of their absolute value is Lebesgue-integrable on G

∞

L (G) set of measurable, essentially bounded functions

scalar product in L2 (G) †

·, · 0,G

norm in L2 (G) †

· 0,G

· 0,p,G (p ∈ [1, ∞]) norm in Lp (G) †

· ∞,G norm in L∞ (G) †

(l ∈ N, p ∈ [1, ∞]) set of l-times weakly di¬erentiable functions

l

Wp (G)

from Lp (G), with derivatives in Lp (G)

· l,p,G (l ∈ N, p ∈ [1, ∞]) norm in Wp (G) †

l

(l ∈ N, p ∈ [1, ∞]) seminorm in Wp (G) †

| · |l,p,G l

A.2. Basic Concepts of Analysis 393

:= W2 (G) (l ∈ N)

H l (G) l

(l ∈ N) scalar product in H l (G) †

·, · l,G

(l ∈ N) norm in H l (G) †

· l,G

(l ∈ N) seminorm in H l (G) †

| · |l,G

·, · 0,h discrete L2 („¦)-scalar product

· 0,h discrete L2 („¦)-norm

L2 (‚G) set of square Lebesgue-integrable functions on the boundary ‚G

1

set of functions from H 1 (G) with vanishing trace on ‚G

H0 (G)

C([0, T ], X) = C 0 ([0, T ], X) set of continuous functions on [0, T ] with values

in the normed space X

C ([0, T ], X) (l ∈ N) set of l-times continuously di¬erentiable functions on

l

[0, T ] with values in the normed space X

Lp ((0, T ), X) (p ∈ [1, ∞]) Lebesgue-space of functions on [0, T ] with values

in the normed space X

†

Convention: In the case G = „¦, this speci¬cation is omitted.

A.2 Basic Concepts of Analysis

A subset G ‚ Rd is called a set of measure zero if, for any number µ > 0, a

countable family of balls Bj with d-dimensional volume µj > 0 exists such

that

∞ ∞

G‚

µj < µ and Bj .

j=1 j=1

Two functions f, g : G ’ R are called equal almost everywhere (in short:

equal a.e., notation: f ≡ g) if the set {x ∈ G : f (x) = g(x)} is of measure

zero.

In particular, a function f : G ’ R is called vanishing almost everywhere

if it is equal to the constant function zero almost everywhere.

A function f : G ’ R is called measurable if there exists a sequence

(fi )i of step functions fi : G ’ R such that fi ’ f for i ’ ∞ almost

everywhere.

In what follows, G denotes a subset of Rd , d ∈ N.

(i) A point x = (x1 , x2 , . . . , xd )T ∈ Rd is called a boundary point of G

if every open neighbourhood (perhaps an open ball) of x contains a

point of G as well as a point of the complementary set R \ G.

(ii) The collection of all boundary points of G is called the boundary of

G and is denoted by ‚G.

(iii) The set G := G ∪ ‚G is called the closure of G.

(iv) The set G is called closed if G = G.

394 A. Appendices

(v) The set G is called open if G © ‚G = ….

(vi) The set G \ ‚G is called the interior of G and is denoted by int G.

A subset G ‚ Rd is called connected if for arbitrary distinct points

x1 , x2 ∈ G there exists a continuous curve in G connecting them.

The set G is called convex if any two points from G can be connected by

a straight-line segment in G.

A nonempty, open, and connected set G ‚ Rd is called a domain in Rd .

By ± = (±1 , . . . , ±d )T ∈ Nd a so-called multi-index is denoted. Multi-

0

indices are a popular tool to abbreviate some elaborate notation. For

example,

d d d

|±| :=

±

±

‚i i

‚ := , ±! := ±i ! , ±i .

i=1 i=1 i=1

The number |±| is called the order (or length) of the multi-index ±.

For a continuous function • : G ’ R, the set supp • := {x ∈ G : •(x) = 0}

denotes the support of •.

A.3 Basic Concepts of Linear Algebra

A square matrix A ∈ Rn,n with entries aij is called symmetric if aij = aji

holds for all i, j ∈ {1, . . . , n}.

A matrix A ∈ Rn,n is called positive de¬nite if x · Ax > 0 for all x ∈

Rn \ {0}.

Given a polynomial p ∈ Pk , k ∈ N0 , of the form

k

with aj ∈ C, j ∈ {0, . . . , k}

aj z j

p(z) =

j=0

and a matrix A ∈ Cn,n , then the following matrix polynomial of A can be

established:

k

aj Aj .

p(A) :=

j=0

Eigenvalues and Eigenvectors

Let A ∈ Cn,n . A number » ∈ C is called an eigenvalue of A if

det(A ’ »I) = 0 .

If » is an eigenvalue of A, then any vector x ∈ Cn \ {0} such that

(” (A ’ »I)x = 0)

Ax = »x

is called an eigenvector of A associated with the eigenvalue ».

A.3. Linear Algebra 395

The polynomial pA (») := det(A ’ »I) is called the characteristic

polynomial of A.

The set of all eigenvalues of a matrix A is called the spectrum of A,

denoted by σ(A).

If all eigenvalues of a matrix A are real, then the numbers »max (A) and

»min (A) denote the largest, respectively smallest, of these eigenvalues.

The number (A) = max»∈σ(A) |»| is called the spectral radius of A.

Norms of Vectors and Matrices

The norm of a vector x ∈ Rn , n ∈ N, is a real-valued function x ’ |x|

satisfying the following three properties:

(i) |x| ≥ 0 for all x ∈ Rn , |x| = 0 ” x = 0 ,

(ii) |±x| = |±| |x| for all ± ∈ R , x ∈ Rn ,

(iii) |x + y| ¤ |x| + |y| for all x, y ∈ Rn .

For example, the most frequently used vector norms are

(a) the maximum norm:

|x|∞ := max |xj | . (A3.1)

j=1...n

p ∈ [1, ∞):

(b) the p -norm,

1/p

n

|x|p := |xj |p . (A3.2)

j=1

The important case p = 2 yields the so-called Euclidean norm:

1/2

n

|x|2 := x2 . (A3.3)

j

j=1

The three most important norms (that is, p = 1, 2, ∞) in Rn are equivalent

in the following sense: The inequalities

√

1

√ |x|2 ¤ |x|∞ ¤ |x|2 ¤ n |x|∞ ,

n

1

|x|1 ¤ |x|∞ ¤ |x|1 ¤ n |x|∞ ,

n

√

1

√ |x|1 ¤ |x|2 ¤ |x|1 ¤ n |x|2

n

are valid for all x ∈ Rn .

The norm of the matrix A ∈ Rn,n is a real-valued function A ’ A

satisfying the following four properties:

A ≥ 0 for all A ∈ Rn,n , A = 0 ” A = 0,

(i)

±A = |±| A for all ± ∈ R , A ∈ Rn,n ,

(ii)

A+B ¤ A + B for all A, B ∈ Rn,n ,

(iii)

AB ¤ A B for all A, B ∈ Rn,n .

(iv)

396 A. Appendices

In comparison with the de¬nition of a vector norm, we include here an

additional property (iv), which is called the submultiplicative property. It

restricts the general set of matrix norms to the practically important class

of submultiplicative norms.

The most common matrix norms are

(a) the total norm:

:= n max |aik | ,

A (A3.4)

G

1¤i,k¤n

(b) the Frobenius norm:

1/2

n

a2

A := , (A3.5)

F ik

i,k=1

(c) the maximum row sum:

n

|aik | ,

A := max (A3.6)

∞

1¤i¤n

k=1

(d) the maximum column sum:

n

|aik | .

A := max (A3.7)

1

1¤k¤n

i=1

All these matrix norms are equivalent. For example, we have

1

¤A ¤A ¤n A p ∈ {1, ∞} ,

A ,

G p G p

n

or

1

¤A ¤A ¤n A

A .

G F G F

n