ńņš. 10 |

2

Ī½:|Ī»Ī½ |ā¤ 1 Ī½:|Ī»Ī½ |ā¤ 1

2 2

2

ļ£« ļ£¶1/2

1ļ£ 1 (k)

c2 ļ£ø

ā¤ = e

Ī½ 2

2 2

Ī½:|Ī»Ī½ |ā¤ 1

2

if the eigenvectors are chosen orthonormal with respect to the Euclidean

scalar product (compare (5.67)). For such an initial error and with ex-

act arithmetic the method would thus have a āsmallā contraction number

independent of the discretization.

For Jacobiā™s method damped by Ļ = 1 this means that if the initial error

2

consists of functions of high frequency only (in the sense of an eigenvector

expansion only of eigenvectors with k or l distant to 1), then the above

considerations hold. But already due to rounding errors we will always

ļ¬nd functions of low frequency in the error such that the above statement

of convergence indeed does not hold, but instead the smoothing property

for the damped Jacobiā™s method is valid: A few steps only lead to a low

reduction of the error but smooth the error in the sense that the parts of

high frequency are reduced considerably.

240 5. Iterative Methods for Systems of Linear Equations

The very idea of the multigrid method lies in the approximative calcula-

tion of this remaining error on a coarse grid. The smooth error can still be

represented on the coarser grid and should be approximated there. Gener-

ally, the dimension of the problem is greatly reduced in this way. Since the

ļ¬nite element discretizations are a central topic of this book, we develop

the idea of multigrid methods for such an example. But it will turn out

that the multigrid method can be used as well for both the ļ¬nite diļ¬erence

and the ļ¬nite volume methods. Multigrid methods have even been success-

fully used in areas other than the discretization of diļ¬erential equations.

Algebraic multigrid methods are generally applicable to systems of linear

equations (5.1) and generate by themselves an abstract analogy of a āgrid

hierarchyā (see, for example, [65]).

5.5.2 Multigrid Method for Finite Element Discretizations

Let Tl = Th be a triangulation that originates from a coarse triangulation

T0 by l applications of a reļ¬nement strategy, for example the strategy of

Section 2.4.1. As we will see, it is not necessary that, for example, in two

space dimensions going from Tk to Tk+1 each triangle will be partitioned

into four triangles. Only the relation

Vk ā‚ Vk+1 , k = 0, . . . , l ā’ 1 ,

has to hold for ļ¬nite-dimensional approximation spaces V0 , V1 , . . . , Vl =

Vh generated by a ļ¬xed ansatz; i.e., the approximation spaces have to be

nested. This holds for all approximation spaces discussed in Section 3.3 if

Tk+1 is still a conforming triangulation and results from Tk by partitioning

of K ā Tk into a possibly varying number of elements of equal kind.

The nodes of Tk , which are the degrees of freedom of the discretization

(possibly multiple in a Hermite ansatz), are denoted by

ak , i = 1, . . . , Mk ,

i

and the corresponding basis functions of Vk are denoted by

Ļ•k , i = 1, . . . , Mk ,

i

with the index k = 0, . . . , l. For a quadratic ansatz on a triangle and Dirich-

let boundary conditions the ak are just the vertices and midpoints of the

i

edges in the interior of the domain. Let the underlying variational equa-

tion (2.21) be deļ¬ned by the bilinear form a and the linear form b on the

function space V . The system of equations to be solved is

Al xl = bl . (5.88)

In addition, we have to consider auxiliary problems

Ak xk = bk

5.5. Multigrid Method 241

for k = 0, . . . , l ā’ 1. For the discretization matrix on each reļ¬nement level

we have, according to (2.34),

(Ak )ij = a Ļ•k , Ļ•k , i, j = 1, . . . , Mk , k = 0, . . . , l ,

j i

and for the right side of the problem to be solved

(bl )i = b Ļ•l , i = 1, . . . , Ml .

i

In Section 2.2, xl is denoted by Ī¾, and bl is denoted by q h .

First we discuss the ļ¬nite element discretization of a variational equa-

tion with symmetric bilinear form, so that in reference to Lemma 2.14 the

Galerkin method to be solved is equivalent to the Ritz method, i.e., to the

minimization of

1

Fl (xl ) := xT Al xl ā’ bT xl .

2l l

Note that l indicates the discretization level and is not an index of a

component or an iteration step.

We distinguish between the function ul ā Vl and the representation

vector xl ā RMl , so that

Ml

xl,i Ļ•l .

ul = (5.89)

i

i=1

For a Lagrange ansatz we have

xl,i = ul al , i = 1, . . . , Ml ,

i

as illustrated by Figure 5.5.

2 1

ui 2

xi = 1

1 1.5

Figure 5.5. ui and xi .

Relation (5.89) deļ¬nes a linear bijective mapping

Pl : RMl ā’ Vl . (5.90)

Thus for z l ā RMl (compare (2.35)),

1T 1

z l Al z l ā’ bT z l = a (Pl z l , Pl z l ) ā’ b(Pl z l ) = F (Pl z l ) ,

Fl (z l ) = l

2 2

where

1

a(u, u) ā’ b(u) for u ā V

F (u) :=

2

242 5. Iterative Methods for Systems of Linear Equations

(k)

Let xl be the kth iterate to the solution of (5.88).

(1) Smoothing step: For ļ¬xed Ī½ ā {1, 2, . . .} calculate

(k+1/2) (k)

= SlĪ½ xl .

xl

Let the corresponding function be:

(k+1/2) (k+1/2)

ā Vl .

ul = Pl xl

(2) Coarse grid correction: Solve (exactly)

(k+1/2)

+ v ā’ min

F ul (5.93)

varying v ā Vlā’1 , with solution vlā’1 . Then set

ĀÆ

= Plā’1 ul + Plā’1 vlā’1 .

(k+1) (k+1/2) (k+1/2)

xl + vlā’1 = xl

ĀÆ ĀÆ

Table 5.5. (k + 1)th step of the two-grid iteration.

is the energy functional for the variational equation.

If xl is an approximation of xl , then the error y l := xl ā’ xl satisļ¬es the

error equation

Al y l = bl ā’ Al xl . (5.91)

This equation is equivalent to the minimization problem

Fl (xl + y l ) = min Fl (xl + y)

yāRMl

and therefore to

F (Pl xl + vl ) = min F (Pl xl + v) , (5.92)

vāVl

with vl = Pl y l .

If the error y l is āsmoothā in the sense that it can be well approxi-

mated also in the lower-dimensional space Vlā’1 , one can solve the error

equation (5.91) approximately as part of an iteration step by solving the

minimization problem (5.92) only on Vlā’1 . The starting condition of a

āsmoothā error will be ensured by the application of a ļ¬xed number of

steps of a smoothing iteration method. Let Sl denote the application of

such a smoothing operation, for example the damped Jacobiā™s method

ā’1

Sl x = x ā’ ĻDl (Al x ā’ bl )

with the diagonal matrix Dl corresponding to Al according to (5.18).

Thus we get the algorithm of the two-grid iteration, whose (k + 1)th step

is described in Table 5.5. Problem (5.93) from Table 5.5 is equivalent to

(compare with Lemma 2.3)

(k+1/2)

for all w ā Vlā’1

a ul + v, w = b(w) (5.94)

5.5. Multigrid Method 243

(1) A priori smoothing: Perform Ī½1 smoothing steps:

(k+1/3) (k)

= SlĪ½1 xl

xl ,

where Ī½1 ā {1, 2, . . .} is ļ¬xed. Let the corresponding function be

(k+1/3) (k+1/3)

ul := Pl xl .

(2) Coarse grid correction: Solve on Vlā’1 the Galerkin discretization

a(ĀÆlā’1 , w) = Ė for all w ā Vlā’1

v b(w) (5.95)

with the bilinear form a and the linear form

(k+1/3)

Ė

b(w) := b(w) ā’ a ul ,w

(a) for l = 1 exactly,

(b) for l > 1 by Āµ steps of a multigrid iteration on level l ā’ 1 for a

and Ė and for the start approximation 0.

b

+ Plā’1 vlā’1 .

(k+2/3) (k+1/3)

Set xl = xl ĀÆ

(3) A posteriori smoothing: Perform Ī½2 smoothing steps

(k+1) (k+2/3)

= SlĪ½2 xl

xl ,

with Ī½2 ā {1, 2, . . .} ļ¬xed.

Table 5.6. (k + 1)th step of the multigrid iteration on level l for bilinear form a

and linear form b.

and thus again to the Galerkin discretization of a variational equation with

Vlā’1 instead of V , with the same bilinear form and with a linear form

deļ¬ned by

(k+1/2)

w ā’ b(w) ā’ a ul for w ā Vlā’1 .

,w

Hence we can ignore the assumption of symmetry for the bilinear form a

and ļ¬nd the approximative solution of the error equation (5.91) on grid

level l ā’ 1 by solving the variational equation (5.94). The equivalent system

of equations will be derived in the following. On the one hand, this prob-

lem has a lower dimension than the original problem, but it also must be

solved for each iteration. This suggests the following recursive procedure:

If we have more than two grid levels, we again approximate this variational

equation by Āµ multigrid iterations; in the same way we treat the hereby

created Galerkin discretization on level l ā’ 2 until level 0 is reached, where

we solve exactly. Furthermore, to conclude each iteration step smoothing

steps should be performed. This leads to the algorithm of the multigrid

iteration. The (k + 1)th step of the multigrid iteration on level l for the

(k)

bilinear form a, linear form b, and starting iteration xl is described in

Table 5.6.

244 5. Iterative Methods for Systems of Linear Equations

In general, Ī½1 = Ī½2 is used. In a convergence analysis it turns out that

only the sum of smoothing steps is important. Despite the recursive deļ¬-

nition of a multigrid iteration we have here a ļ¬nite method, because the

level 0 is reached after at most l recursions, where the auxiliary problem

will be solved exactly. For Āµ usually only the values Āµ = 1 or Āµ = 2 are

used. The terms V-cycle for Āµ = 1 and W-cycle for Āµ = 2 are commonly

used, because for an iteration, the sequence of levels on which operations

are executed have the shape of these letters (see Figure 5.6).

for l = 2 :

Level

Āµ=1 Āµ=2

2 o o o o

1 o o o o o

0 o o o

for l = 3 :

Level

3 o o o o

2 o o o o o

1 o o o o o o o o

0 o o o o o

Figure 5.6. Grid levels for the V-cycle (Āµ = 1) and the W-cycle (Āµ = 2).

The problems in (5.94) and (5.95) (see Table 5.6) have the form

for all w ā Vlā’1 ,

a (u + v, w) = b(w) (5.96)

where v ā Vlā’1 is unknown and u ā Vl is known. An equivalent system of

equations arises by inserting the basis functions Ļ•lā’1 , j = 1, . . . , Mlā’1 ,

j

for w and an appropriate representation for v. If we again take the

representation with respect to Ļ•lā’1 , we get as in (2.34)

j

ā’1

Alā’1 Plā’1 v = dlā’1 . (5.97)

Here the residual dk ā RMk of u on the diļ¬erent levels k = 0, . . . , l is

deļ¬ned by

dk,i := b Ļ•k ā’ a u, Ļ•k , i = 1, . . . , Mk .

i i

5.5. Multigrid Method 245

We now develop an alternative representation for (5.97) and the coarse

grid correction for possible generalizations beyond the Galerkin approxima-

tions. Therefore, let R ā RMlā’1 ,Ml be the matrix that arises through the

unique representation of the basis functions Ļ•lā’1 with respect to the basis

j

l

Ļ•i , which means the elements rji of R are determined by the equations

Ml

Ļ•lā’1 rji Ļ•l ,

= j = 1, . . . , Mlā’1 .

i

j

i=1

Then (5.96) is equivalent to

a(v, w) = b(w) ā’ a(u, w) for all w ā Vlā’1

ļ£« ļ£¶

Mlā’1

ā” aļ£ Ļ•lā’1 , Ļ•lā’1 ļ£ø = b Ļ•lā’1 ā’ a u, Ļ•lā’1 , j = 1, . . . , Mlā’1

ā’1

Plā’1 v s j j j

s

s=1

Mlā’1 Ml Ml Ml

ā’1

ā” rji b Ļ•l ā’ a u, Ļ•l

rst Ļ•l rji Ļ•l

Plā’1 v s a , =

t i i i

s=1 t=1 i=1 i=1

Mlā’1 Ml

ā’1

ā” rji a Ļ•l , Ļ•l rst Plā’1 v = (Rdl )j , j = 1, . . . , Mlā’1 .

t i s

s=1 i,t=1

Hence the system of equations has the form

ā’1

RAl RT Plā’1 v = Rdl . (5.98)

The matrix R is easy to calculate for a node-based basis Ļ•l satisfying

i

Ļ•i aj = Ī“ij , since in this case we have for v ā Vl ,

l l

Ml

v al Ļ•l ,

v= i i

i=1

and therefore in particular,

Ml

Ļ•lā’1 Ļ•lā’1 al Ļ•l

= i i

j j

i=1

and thus

rji = Ļ•lā’1 al .

i

j

For the linear ansatz in one space dimension with Dirichlet boundary

1

conditions (i.e., with V = H0 (a, b) as basic space) this means that

ļ£«1 ļ£¶

11

2 2

ļ£¬ ļ£·

1

11

ļ£¬ ļ£·

R=ļ£¬ ļ£·.

2 2

(5.99)

ļ£¬ ļ£·

..

ļ£ ļ£ø

.

1 1

1

2 2

246 5. Iterative Methods for Systems of Linear Equations

The representation (5.98) can also be interpreted in this way:

Due to Vlā’1 ā‚ Vl the identify deļ¬nes a natural prolongation from Vlā’1 to

Vl , which means that

p : Vlā’1 ā’ Vl , v ā’ v ,

Ė

as illustrated by Figure 5.7.

1 1

1 1

2 2

Figure 5.7. Prolongation.

This prolongation corresponds to a prolongation p from RMlā’1 to RMl ,

the canonical prolongation, through the transition to the representation

vectors (5.90). It is given by

p := Plā’1 Plā’1 , (5.100)

since for xlā’1 ā RMlā’1 , p can be composed as follows:

p

Ė

xlā’1 ā’ Plā’1 xlā’1 ā’ Plā’1 xlā’1 ā’ Plā’1 Plā’1 xlā’1 .

Obviously, p is linear and can be identiļ¬ed with its matrix representation

in RMl ,Mlā’1 . Then

p = RT (5.101)

holds, because

Mlā’1 Ml Mlā’1

yj Ļ•lā’1 yj rji Ļ•l ,

Plā’1 y = = i

j

j=1 i=1 j=1

i.e., RT y = Plā’1 (Plā’1 y) for any y ā RMlā’1 .

(l)

In the following RMl will be endowed with a scalar product Ā·, Ā· , which

is an Euclidean scalar product scaled by a factor Sl ,

Ml

(l)

xl , y l := Sl xl,i yl,i . (5.102)

i=1

The scaling factor is to be chosen such that for the induced norm Ā· and

l

the L2 (ā„¦)-norm on Vl ,

ā¤ xl ā¤ C2 Pl xl

C1 Pl xl (5.103)

0 l 0

for x ā RMl , l = 0, 1, . . ., with constants C1 , C2 independent of l: If the

triangulations are members of a regular and quasi-uniform family Th (see

5.5. Multigrid Method 247

Deļ¬nition 3.28), then in d space dimensions one can choose Sl = hd , with

l

hl being the maximal diameter of K ā Tl (see Theorem 3.43).

Let r : RMl ā’ RMlā’1 be deļ¬ned by

r = pā— , (5.104)

with the adjoint pā— deļ¬ned with respect to the scalar products Ā·, Ā· (lā’1)

and Ā·, Ā· (l) ; that is,

(lā’1) (lā’1) (l)

= pā— xl , y lā’1

r xl , y lā’1 = xl , p y lā’1 .

If p is the canonical prolongation, then r is called the canonical restriction.

For the representation matrices,

Slā’1

r = pT = R . (5.105)

Sl

In example (5.102) for d = 2 with hl = hlā’1 /2 we have Slā’1 /Sl = 1/4. Due

to Pl p = Plā’1 , the canonical restriction of RMl on RMlā’1 satisļ¬es

rRl = Rlā’1 ,

where Rl : Vl ā’ RMl is deļ¬ned as the adjoint of Pl ,

(l)

for all xl ā RMl , vl ā Vl ,

Pl xl , vl = xl , Rl vl

0

because for any y lā’1 ā RMlā’1 and for vlā’1 ā Vlā’1 ā‚ Vl ,

(lā’1) (l)

rRl vlā’1 , y lā’1 = Rl vlā’1 , pylā’1 = vlā’1 , Pl pylā’1 0

(lā’1)

= vlā’1 , Plā’1 y lā’1 = Rlā’1 vlā’1 , y lā’1 .

0

Using (5.105) we see the equivalence of equation (5.98) to

(rAl p)y lā’1 = rdl . (5.106)

Ė Ė

Setting v := Plā’1 y lā’1 for a perhaps only approximative solution y lā’1 of

ā’1

(5.106), the coarse grid correction will be ļ¬nished by addition of Pl v. Due

to

Plā’1 v = Plā’1 Plā’1 Plā’1 v = p Plā’1 v ,

ā’1 ā’1

the coarse grid correction is

(k+2/3) (k+1/3)

xl = xl + p(Ė lā’1 ) .

y

The above-mentioned facts suggest the following structure of a general

multigrid method: For discretizations deļ¬ning a hierarchy of discrete

problems,

Al xl = bl ,

one needs prolongations

p : RMkā’1 ā’ RMk

248 5. Iterative Methods for Systems of Linear Equations

and restrictions

r : RMk ā’ RMkā’1

Ė

for k = 1, . . . , l and the matrices Akā’1 for the error equations. The coarse

grid correction steps (5.93) and (5.95) hence take the following form:

Solve (with Āµ steps of the multigrid method)

(k+1/3)

Alā’1 y lā’1 = r bl ā’ Al xl

Ė

and set

(k+2/3) (k+1/3)

xl = xl + py lā’1 .

The above choice

Ė

Alā’1 = rAl p

is called the Galerkin product. For Galerkin approximations this coincides

with the discretization matrix of the same type on the grid of level l ā’ 1 due

to (5.97). This is also a common choice for other discretizations and then

an alternative to the Galerkin product. In view of the choice of p and r we

should observe the validity of (5.104). An interpolational deļ¬nition of the

prolongation on the basis of (ļ¬nite element) basis functions as for example

(5.101) (see also example (5.99)) is also common in other discretizations. In

more diļ¬cult problems, as for example those with (dominant) convection

in addition to diļ¬usive transport processes, nonsymmetric problems arise

with a small constant of V -ellipticity. Here the use of matrix-dependent,

that means Al -dependent, prolongations and restrictions is recommended.

5.5.3 Eļ¬ort and Convergence Behaviour

In order to judge the eļ¬ciency of a multigrid method the number of opera-

tions per iteration and the number of iterations (required to reach an error

level Īµ, see (5.4)) has to be estimated. Due to the recursive structure, the

ļ¬rst number is not immediately clear. The aim is to have only the optimal

amount of O(Ml ) operations for sparse matrices. For this the dimensions

of the auxiliary problems have to decrease suļ¬ciently. This is expressed by

the following:

There exists a constant C > 1 such that

Mlā’1 ā¤ Ml /C for l ā N . (5.107)

Hence we assume an inļ¬nite hierarchy of problems and/or grids, which

also corresponds to the asymptotic point of view of a discretization from

Section 3.4. Relation (5.107) is thus a condition for a reļ¬nement strategy.

For the model problem of the Friedrichsā“Keller triangulation of a rectangle

(see Figure 2.9) in the case of a regular āredā reļ¬nement we have hl =

hlā’1 /2. Thus C = 4, and for analogous constructions in d space dimensions

5.5. Multigrid Method 249

C = 2d . The matrices that appear should be sparse, so that for level l the

following holds:

smoothing step = CS Ml operations,

error calculation and restrictions = CD Ml operations,

prolongation and correction = CC Ml operations.

Then we can prove the following (see [16, p. 326]):

If the number Āµ of multigrid steps in the recursion satisļ¬es

Āµ<C, (5.108)

then the number of operations for an iteration step for a problem on level

l can be estimated by

C(Ī½)Ml . (5.109)

Here Ī½ is the number of a priori and a posteriori smoothing steps and

Ī½CS + CD + CS

+ O (Āµ/C)l .

C(Ī½) =

1 ā’ Āµ/C

The requirement (5.108) will be satisļ¬ed in general through the restriction

to Āµ = 1, Āµ = 2. Analogously, the memory requirement is O(Ml ), since

l

C

Mk ā¤ Ml .

C ā’1

k=0

Whether this larger eļ¬ort (of equal complexity) in comparison to other

methods discussed is justiļ¬ed will be decided by the rate of convergence.

The multigrid method is a linear stationary method. The iteration matrix

MlT GM of the two-grid method results from

(k+1/2) (k)

= SlĪ½ xl

xl ,

+ p Aā’1 r bl ā’ Al xl

(k+1) (k+1/2) (k+1/2)

xl = xl lā’1

to

MlT GM = (I ā’ pAā’1 rAl )SlĪ½ . (5.110)

lā’1

Also, the consistency of the method follows immediately if the smoothing

iteration is consistent.

The analysis of the convergence of the multigrid method can be reduced

to the analysis of the two-grid method, since the iteration matrix is a

modiļ¬cation of MlT GM (see [16, p. 328]). For a large class of a priori and a

posteriori smoothing operators as well as of restrictions and prolongations

it can be shown (see [16, p. 347]) that there exists a constant ā (0, 1)

independent of the discretization parameter hl such that (MT GM ) ā¤ .

Combined with (5.109) this shows that multigrid methods are optimal in

their complexity. This also shows their potential superiority compared with

all other methods described.

250 5. Iterative Methods for Systems of Linear Equations

In the following we will only indicate the schematic procedure to prove

this assertion. It is suļ¬cient to prove the following two properties, where

the spectral norm is used as the matrix norm, that is, the matrix norm

that is induced by the Euclidean vector norm.

(1) Smoothing property:

CS

Al SlĪ½ ā¤

There exists CS > 0 : Al .

Ī½

(2) Approximation property:

Aā’1 ā’ pAā’1 r ā¤ CA Al ā’1

There exists CA > 0 : . (5.111)

l lā’1

Due to

MT GM = Aā’1 ā’ pAā’1 r Al SlĪ½ ,

l lā’1

we can conclude that

CS CA

MT GM ā¤ Aā’1 ā’ pAā’1 r Al SlĪ½ ā¤ ,

l lā’1

Ī½

which means that for suļ¬ciently large Ī½,

MT GM ā¤ <1

with independent of l.

The smoothing property is of an algebraic nature, but for the proof of

the approximation property we will use ā” at least indirectly ā” the original

variational formulation of the boundary value problem and the correspond-

ing error estimate. Therefore, we discuss only the smoothing property for,

as an example, the relaxed Richardson method for a symmetric positive

deļ¬nite matrix Al , i.e.,

1

Sl = Il ā’ ĻAl Ļā

with 0, .

Ī»max (Al )

Let {z i }Ml be an orthonormal basis of eigenvectors of Al . For any initial

i=1

Ml

vector x(0) represented in this basis as x(0) = i=1 ci z i it follows that

(compare (5.68))

Ml Ml

2 ā’2

ā’ (Ī»i Ļ)2 (1 ā’ Ī»i Ļ)2Ī½ c2

Al SlĪ½ x(0) Ī»2 (1 Ī»i Ļ)2Ī½ c2

= =Ļ

i i i

i=1 i=1

2 Ml

ā’2

ā¤Ļ max Ī¾(1 ā’ Ī¾) Ī½

c2 .

i

Ī¾ā[0,1]

i=1

The function Ī¾ ā’ Ī¾(1 ā’ Ī¾)Ī½ has its maximum at Ī¾max = (Ī½ + 1)ā’1 ; thus

Ī½ Ī½+1

1 1 1 Ī½ 1

Ī¾max (1 ā’ Ī¾max ) = 1ā’ ā¤

Ī½

= .

Ī½ +1 Ī½ +1 Ī½ Ī½+1 eĪ½

5.6. Nested Iterations 251

Hence

1

Al SlĪ½ x(0) ā¤ x(0) ,

ĻeĪ½

which implies

1

Al SlĪ½ ā¤ .

ĻeĪ½

Since the inclusion Ļ ā (0, 1/Ī»max(Al )] can be written in the form Ļ =

Ļ/ Al with Ļ ā (0, 1], we have CS = 1/(Ļe).

The approximation property can be motivated in the following way. The

ļ¬ne grid solution xl of Al xl = dl is replaced in the coarse grid correction

by pxlā’1 from Alā’1 xlā’1 = dlā’1 := rdl . Therefore, pxlā’1 ā Aā’1 dl should

l

hold. The formulation (5.111) thus is just a quantitative version of this

requirement. Since in the symmetric case Al ā’1 is simply the reciprocal

value of the largest eigenvalue, (3.140) in Theorem 3.45 establishes the

relation to the statements of convergence in Section 3.4. For a more exact

analysis of convergence and a more extensive description of this topic we

refer to the cited literature (see also [17]).

Exercises

5.12 Determine the prolongation and restriction according to (5.101) and

(5.104) for the case of a linear ansatz on a Friedrichsā“Keller triangulation.

5.13 Prove the consistency of the two-grid method (5.110) in the case of

the consistent smoothing property.

5.6 Nested Iterations

As in Section 5.5 we assume that besides the system of equations

Al xl = bl

with Ml unknowns, there are given analogous low-dimensional systems of

equations

k = 0, . . . , l ā’ 1 ,

Ak xk = bk , (5.112)

with Mk unknowns, where M0 < M1 < Ā· Ā· Ā· < Ml . Let all systems of

equations be an approximation of the same continuous problem such that

an error estimate of the type

u ā’ Pl xl ā¤ CA hĪ±

l

holds, with Pl according to (5.90) and Ī± > 0. Here Ā· is a norm on the

basic space V , and the constant CA generally depends on the solution u

252 5. Iterative Methods for Systems of Linear Equations

of the continuous problem. The discretization parameter hl determines the

dimension Ml : In the simplest case of a uniform reļ¬nement, hd ā¼ 1/Ml

l

holds in d space dimensions. One may also expect that for the discrete

solution,

pxkā’1 ā’ xk ā¤ C1 CA hĪ± , k = 1, . . . , l ,

k k

holds with a constant C1 > 0. Here Ā· k is a norm on RMk , and the

mapping p = pkā’1,k : RMkā’1 ā’ RMk is a prolongation, for example the

canonical prolongation introduced in Section 5.5. In this case the estimate

can be rigorously proven with the deļ¬nition of the canonical prolongation

ā’1

p = Pk Pkā’1 :

ā’1

pxkā’1 ā’ xk Pk Pkā’1 xkā’1 ā’ Pk xk

=

k k

ā’1

ā¤ Pkā’1 xkā’1 ā’ Pk xk

Pk L[V ,RMk ]

k

ā’1

ā¤ ā¤ C1 CA hĪ±

CA hĪ± + CA hĪ±

Pk L[V ,RMk ] k kā’1 k

k

with

Ī±

hjā’1

Pjā’1 L[Vj ,RMj ]

C1 = max 1+ .

hj

j=1,...,l

Let the system of equations be solved with an iterative method given by

the ļ¬xed-point mapping Ī¦k , k = 0, . . . , l, which means that xk according

to (5.112) satisļ¬es xk = Ī¦k (xk , bk ). Then it is suļ¬cient to determine an

Ė

iterate xl with an accuracy

xl ā’ xl l ā¤ CA hĪ±

Ė

Ė (5.113)

l

Ė

with CA := CA / Pl L[RMl ,V ] , because then we also have

Pl xl ā’ Pl xl ā¤ CA hĪ± .

Ė l

If one does not have a good initial iterate from the concrete context, the

algorithm of nested iterations explained in Table 5.7 can be used. It is

indeed a ļ¬nite process.

The question is how to choose the iteration numbers mk such that (5.113)

ļ¬nally holds, and whether the arising overall eļ¬ort is acceptable. An answer

to this question is provided by the following theorem:

Theorem 5.19 Let the iterative method Ī¦k have the contraction number

k with respect to Ā· k . Assume that there exist constants C2 , C3 > 0 such

that

ā¤ C2 ,

p L[RMkā’1 ,RMk ]

ā¤ C3 hk ,

hkā’1

for all k = 1, . . . , l. If the iteration numbers mk for the nested iterations

are chosen in such a way that

ā¤ 1/(C2 C3 + C1 Pl ) ,

mk Ī±

(5.114)

k

5.6. Nested Iterations 253

Choose mk , k = 1, . . . , l.

Ė

Let x0 be an approximation of x0 ,

for example x0 = x0 = Aā’1 b0 .

Ė 0

For k = 1, . . . , l:

(0)

Ė Ė

xk := p xkā’1 .

Perform mk iterations:

(i) (iā’1)

Ė Ė

xk := Ī¦k xk , bk , i = 1, . . . , mk .

(mk )

Ė Ė

Set xk := xk .

Table 5.7. Nested Iteration.

then

xk ā’ xk ā¤ CA hĪ± ,

Ėk

Ė k

for all k = 1, . . . , l, provided that this estimate holds for k = 0.

Proof: The proof is given by induction on k. Assume that the assertion

is true for k ā’ 1. This induces

xk ā’ xk ā¤ pĖ kā’1 ā’ xk

mk

Ė x

k k k

ā¤ ( p(Ė kā’1 ā’ xkā’1 ) + pxkā’1 ā’ xk

mk

x k)

k

k

ā¤ Ė

mk

C2 CA hĪ± + C1 CA hĪ±

kā’1 k

k

ā¤ Ėk

mk

(C2 C3 + C1 Pl ) CA hĪ± .

Ī±

k

2

Theorem 5.19 allows the calculation of the necessary number of iterations

ā’1

for the inner iteration from the norms p L[RMkā’1 ,RMk ] , Pk L[Vk ,RMk ] and

the constants hhk for k = 1, . . . , l, as well as the order of convergence Ī±

kā’1

of the discretization.

In order to estimate the necessary eļ¬ort according to (5.114) more ex-

actly, the dependence of k of k must be known. In the following we consider

only the situation, known as the multigrid method, of a method of optimal

complexity

ā¤ < 1.

k

Here, in contrast to other methods, the number of iterations can be cho-

sen constant (mk = m for all k = 1, . . . , l). If, furthermore, the estimate

(5.107) holds with the constant C, then analogously to the consideration

in Section 5.5 the total number of operations for the nested iteration can

254 5. Iterative Methods for Systems of Linear Equations

be estimated by

C

m CMl .

C ā’1

Here CMk is the number of operations for an iteration with the iteration

method Ī¦k .

In the model problem of the Friedrichsā“Keller triangulation with uniform

reļ¬nement we have C/(C ā’ 1) = 4/3 and C3 = 2. For Ā· = Ā· 0 as basic

norm, Ī± = 2 is a typical case according to Theorem 3.37. The existence of

the constant C2 will hereby ļ¬nally be ensured consistently by the condition

(5.103), observing (5.100). Assuming also that the constants C1 , C2 , Pl

are āsmallā and the iteration method has a āsmallā contraction number ,

only a small number of iterations m is necessary, in the ideal case m = 1.

At least in this situation we can count on only a small increase of the

necessary eļ¬ort through the process of nested iterations, which provides an

Ė

āappropriateā approximation xk on all levels k of discretization.

Finally, it is to be observed that the sequence of the discrete problems

has to be deļ¬ned only during the process of the nested iteration. This oļ¬ers

the possibility to combine it with a posteriori error estimators as discussed

in Section 4.2, in order to develop a grid Tk+1 on which the discrete problem

of level k + 1 is determined, on the basis of xk as a reļ¬nement of Tk .

Ė

6

The Finite Volume Method

Finite volume methods are widely applied when diļ¬erential equations in

divergence form (cf. Section 0.5) or diļ¬erential equations involving such

diļ¬erential expressions (for example, parabolic diļ¬erential equations) are to

be solved numerically. In the class of second-order linear elliptic diļ¬erential

equations, expressions of the form

Lu := ā’ā Ā· (K āu ā’ c u) + r u = f (6.1)

are typical (cf. (0.33)), where

K : ā„¦ ā’ Rd,d , c : ā„¦ ā’ Rd , r, f : ā„¦ ā’ R .

The corresponding āparabolic versionā is

ā‚u

+ Lu = f

ā‚t

and will be treated in Chapter 7.

First-order partial diļ¬erential equations such as the classical conservation

laws

ā Ā· q(u) = 0 ,

where q : R ā’ Rd is a nonlinear vector ļ¬eld depending on u, or higher-order

partial diļ¬erential equations (such as the biharmonic equation (3.36)), or

even systems of partial diļ¬erential equations can be successfully discretized

by the ļ¬nite volume method.

In correspondence to the comparatively large class of problems that can

be treated by the ļ¬nite volume method, there are rather diļ¬erent sources

256 6. Finite Volume Method

1960 Forsythe and Wasow computation of neutron diļ¬usion

1961 MarĖuk

c computation of nuclear reactors

1971 McDonald ļ¬‚uid mechanics

1972 MacCormack and Paullay ļ¬‚uid mechanics

1973 Rizzi and Inouye ļ¬‚uid mechanics in 3D

1977 Samarski integro-interpolation method,

balance method

.

.

.

1979 Jameson ļ¬nite volume method

1984 Heinrich integro-balance method,

generalized ļ¬nite diļ¬erence method

.

.

.

1987 Bank and Rose box method

.

.

.

Table 6.1. Some sources of the ļ¬nite volume method.

originating mainly from practical applications. Some of these sources are

listed in Table 6.1. In contrast to ļ¬nite diļ¬erence or ļ¬nite element methods,

the theoretical understanding of the ļ¬nite volume method remained at an

early stage for a long time; only in recent years has essential progress been

noted.

The ļ¬nite volume method can be viewed as a discretization method of

its own right. It includes ideas from both ļ¬nite diļ¬erence and ļ¬nite element

methods. So in the literature approaches can be found that interpret it as

a āgeneralized ļ¬nite diļ¬erence methodā or rather as a variant of the ļ¬nite

element method. In this chapter, we will consider only equations of the

type (6.1).

6.1 The Basic Idea of the Finite Volume Method

Now we will describe the fundamental steps in the derivation of the ļ¬nite

volume method. For simplicity, we restrict ourselves to the case d = 2 and

r = 0. Furthermore, we set q(u) := ā’K āu + c u. Then equation (6.1)

becomes

ā Ā· q(u) = f . (6.2)

In order to obtain a ļ¬nite volume discretization, the domain ā„¦ will be

subdivided into M subdomains ā„¦i such that the collection of all those

subdomains forms a partition of ā„¦, that is:

(1) each ā„¦i is an open, simply connected, and polygonally bounded set

without slits,

6.1. Basics 257

(2) ā„¦i ā© ā„¦j = ā… (i = j),

(3) āŖM ā„¦i = ā„¦ .

i=1

These subdomains ā„¦i are called control volumes or control domains.

Without going into more detail we mention that there also exist ļ¬nite

volume methods with a well-deļ¬ned overlapping of the control volumes

(that is, condition 2 is violated).

The next step, which is in common with all ļ¬nite volume methods, con-

sists in integrating equation (6.2) over each control volume ā„¦i . After that,

Gaussā™s divergence theorem is applied:

Ī½ Ā· q(u) dĻ = i ā {1, . . . , M } ,

f dx ,

ā‚ā„¦i ā„¦i

where Ī½ denotes the outer unit normal to ā‚ā„¦i . By the ļ¬rst condition of the

partition, the boundary ā‚ā„¦i is formed by straight-line segments Ī“ij (j =

1, . . . , ni ), along which the normal Ī½|Ī“ij =: Ī½ij is constant (see Figure 6.1).

So the line integral can be decomposed into a sum of line integrals from

which the following equation results:

ni

Ī½ij Ā· q(u) dĻ = f dx . (6.3)

Ī“ij ā„¦i

j=1

Ī½i 2 Ī½i 1

ā„¦i

Ī½i 3

Ī½i 5

Ī½i 4

Figure 6.1. A control volume.

Now the integrals occurring in (6.3) have to be approximated. This can

be done in very diļ¬erent ways, and so diļ¬erent ļ¬nal discretizations are

obtained.

In general, ļ¬nite volume methods can be distinguished by the following

criteria:

(1) the geometric shape of the control volumes ā„¦i ,

(2) the position of the unknowns (āproblem variablesā) with respect to

the control volumes,

258 6. Finite Volume Method

(3) the approximation of the boundary (line (d = 2) or surface (d = 3))

integrals.

Especially the second criterion divides the ļ¬nite volume methods into two

large classes: the cell-centred and the cell-vertex ļ¬nite volume methods. In

the cell-centred methods, the unknowns are associated with the control

volumes (for example, any control volume corresponds to a function value

at some interior point (e.g., at the barycentre)). In the cell-vertex methods,

the unknowns are located at the vertices of the control volumes. Sometimes,

instead of the ļ¬rst-mentioned class a subdivision into two classes, the so-

called cell-centred and node-centred methods, is considered. The diļ¬erence

is whether the problem variables are assigned to the control volumes or,

given the problem variables, associated control volumes are deļ¬ned.

Example 6.1 Consider the homogeneous Dirichlet problem for the Pois-

son equation on the unit square:

ā’āu = in ā„¦ = (0, 1)2 ,

f

u= 0 on ā‚ā„¦ .

Problem variables:

aj2

Function values at the nodes ai

of a square grid with mesh width

aj3 ai aj1

h>0

aj4 Control volumes:

ā„¦i := {x ā ā„¦ : |x ā’ ai |ā < h }

2

Figure 6.2. Problem variables and control volumes in a cell-centred ļ¬nite volume

method.

For an inner control volume ā„¦i (i.e., ai ā ā„¦), equation (6.3) takes the

form

4

ā’ Ī½ijk Ā· āu dĻ = f dx ,

Ī“ijk ā„¦i

k=1

where Ī“ijk := ā‚ā„¦i ā© ā‚ā„¦jk . A closer look at the directional derivatives shows

that

Ī½ij1 Ā· āu = ā‚1 u , Ī½ij2 Ā· āu = ā‚2 u ,

Ī½ij3 Ā· āu = ā’ā‚1 u , Ī½ij4 Ā· āu = ā’ā‚2 u .

i.e. they are just partial derivatives with respect to the ļ¬rst or the second

variable on the corresponding parts of the boundary.

6.1. Basics 259

Approximating the integrals on Ī“ijk by means of the midpoint rule and

replacing the derivatives by diļ¬erence quotients, we have

4 4

ai + ajk

ā’ Ī½ijk Ā· āu dĻ ā ā’ Ī½ijk Ā· āu h

2

Ī“ijk

k=1 k=1

u(aj1 )ā’u(ai ) u(aj2 )ā’u(ai ) u(ai )ā’u(aj3 ) u(ai )ā’u(aj4 )

āā’ ā’ ā’

+ h

h h h h

4

= 4 u(ai ) ā’ u(ajk ) .

k=1

Thus, we obtain exactly the expression that results from the application of

a ļ¬nite element method with continuous, piecewise linear ansatz and test

functions on a Friedrichsā“Keller triangulation (cf. Figure 2.9).

Furthermore, if we approximate the integral ā„¦i f dx by f (ai )h2 , we see

that this term coincides with the trapeziodal rule applied to the right-hand

side of the mentioned ļ¬nite element formulation (cf. Lemma 2.13).

Actually, it is no accident that both discretization methods lead to the

same algebraic system. Later on we will prove a more general result to

conļ¬rm the above observation.

The boundary control volumes are treated as follows:

If ai ā ā‚ā„¦, then parts of the boundary ā‚ā„¦i lie on ā‚ā„¦. At these nodes,

the Dirichlet boundary conditions already prescribe values of the unknown

function, and so there is no need to include the boundary control volumes

into the balance equations (6.3).

A detailed description for the case of ļ¬‚ux boundary conditions will be

given later, in Section 6.2.4; see (6.23).

Example 6.2 We consider the same boundary value problem as in

Example 6.1.

Problem variables:

Function values at the nodes ai

of a square grid with mesh width

h>0

ā„¦i

Control volumes:

Subsquares of the grid

Figure 6.3. Problem variables and control volumes in a cell-vertex ļ¬nite volume

method.

260 6. Finite Volume Method

In the interior of ā„¦, the resulting discretization yields a 12-point stencil

(in the terminology of ļ¬nite diļ¬erence methods).

Remark 6.3 In the ļ¬nite volume discretization of systems of partial dif-

ferential equations (resulting from ļ¬‚uid mechanics, for example), both

methods are used simultaneously for diļ¬erent variables; see Figure 6.4.

. . . ..

O O O O O

. . . . : problem variable of type 1

O O O O O

. . . . : problem variable of type 2

O

O O O O O

. . . .

O O O O O

O O O O O

Figure 6.4. Finite volume discretization of systems of partial diļ¬erential

equations.

Assets and Drawbacks of the Finite Volume Method

Assets:

ā¢ Flexibility with respect to the geometry of the domain ā„¦ (as in ļ¬nite

element methods).

ā¢ Admissibility of unstructured grids (as in ļ¬nite element methods,

important for adaptive methods).

ā¢ Simple assembling.

ā¢ Conservation of certain laws valid for the continuous problem

(for example, conservation laws or maximum principles). This

property is important in the numerical solution of diļ¬erential equa-

tions with discontinuous coeļ¬cients or of convection-dominated

diļ¬usion-convection equations (see Section 6.2.4).

ā¢ Easy linearization of nonlinear problems (simpler than in ļ¬nite

element methods (Newtonā™s method)).

ā¢ Simple discretization of boundary conditions (as in ļ¬nite element

methods, especially a ānaturalā treatment of Neumann or mixed

boundary conditions).

ā¢ In principle, no restriction of the spatial dimension d of the domain

ā„¦.

6.1. Basics 261

Drawbacks:

ā¢ Smaller ļ¬eld of applications in comparison with ļ¬nite element or ļ¬nite

diļ¬erence methods.

ā¢ Diļ¬culties in the design of higher order methods (no so-called p-

version available as in the ļ¬nite element method).

ā¢ In higher spatial dimensions (d ā„ 3), the construction of some classes

or types of control volumes may be a complex task and thus may lead

to a time-consuming assembling.

ā¢ Diļ¬cult mathematical analysis (stability, convergence, . . . ).

Exercises

6.1 Given the boundary value problem

ā’(au ) = 0 in (0, 1) , u(0) = 1 , u(1) = 0 ,

with piecewise constant coeļ¬cients

ĪŗĪ± , x ā (0, Ī¾) ,

a(x) :=

x ā (Ī¾, 1) ,

Ī±,

where Ī±, Īŗ are positive constants and Ī¾ ā (0, 1) \ Q :

(a) What is the weak solution u ā H 1 (0, 1) of this problem?

(b) For general āsmoothā coeļ¬cients a, the diļ¬erential equation is

obviously equivalent to

ā’au ā’ a u = 0 .

Therefore, the following discretization is suggested:

uiā’1 ā’ 2ui + ui+1 ai+1 ā’ aiā’1 ui+1 ā’ uiā’1

ā’ai ā’ = 0,

2

h 2h 2h

where an equidistant grid with the nodes xi = ih (i = 0, . . . , N + 1)

and ai := a(xi ), ui :ā u(xi ) is used.

This discretization is also formally correct in the given situation of

discontinuous coeļ¬cients. Find the discrete solution (ui )N in this

i=1

case.

(c) Under what conditions do the values ui converge to u(xi ) for h ā’ 0?

262 6. Finite Volume Method

6.2 The Finite Volume Method for Linear Elliptic

Diļ¬erential Equations of Second Order on

Triangular Grids

In this section we will explain the development and the analysis of a ļ¬nite

volume method of ācell-centredā type for a model problem. Here, ā„¦ ā‚ R2

is a bounded, simply connected domain with a polygonal boundary, but

without slits.

6.2.1 Admissible Control Volumes

The Voronoi Diagram

By {ai }iāĪ ā‚ ā„¦ we denote a consecutively numbered point set that includes

all vertices of ā„¦, where Ī is the corresponding set of indices. Typically, the

points ai are placed at those positions where the values u(ai ) of the exact

solution u are to be approximated. The convex set

ā„¦i := x ā R2 |x ā’ ai | < |x ā’ aj |

Ė for all j = i

is called the Voronoi polygon (or Dirichlet domain, Wignerā“Seitz cell,

Ė

Thiessen polygon, . . . ). The family ā„¦i iāĪ is called the Voronoi diagram

of the point set {ai }iāĪ .

. . .

.. boundary of ā„¦

. ā¼

.. boundary of ā„¦i

Figure 6.5. Voronoi diagram.

The Voronoi polygons are convex, but not necessarily bounded, sets (con-

sider the situation near the boundary in Figure 6.5). Their boundaries are

polygons. The vertices of these polygons are called Voronoi vertices.

It can be shown that at any Voronoi vertex at least three Voronoi poly-

gons meet. According to this property, Voronoi vertices are classiļ¬ed into

regular and degenerate Voronoi vertices: In a regular Voronoi vertex, the

boundaries of exactly three Voronoi polygons meet, whereas a degenerate

Voronoi vertex is shared by at least four Voronoi polygons. In the latter

case, all the corresponding nodes ai are located at some circle (they are

ācocyclicā, cf. Figure 6.6).

6.2. Finite Volume Method on Triangular Grids 263

.a2

. . (a1 - a4 are cocyclic)

a3 a5

.

. a1

a4

Figure 6.6. Degenerate and regular Voronoi vertex.

Now the elements ā„¦i (control volumes) of the partition of ā„¦ required for

the deļ¬nition of the ļ¬nite volume method can be introduced as follows:

ā„¦i := ā„¦i ā© ā„¦ , i ā Ī.

Ė

As a consequence, the domains ā„¦i need not necessarily be convex if ā„¦ is

nonconvex (cf. Figure 6.5).

Furthermore, the following notation will be used:

j ā Ī \ {i} : ā‚ā„¦i ā© ā‚ā„¦j = ā… , i ā Ī,

Īi :=

for the set of indices of neighbouring nodes,

ā‚ā„¦i ā© ā‚ā„¦j , j ā Īi , for a joint piece of the

Ī“ij :=

boundaries of neighbouring control volumes,

mij for the length of Ī“ij .

The dual graph of the Voronoi diagram is deļ¬ned as follows:

Any pair of points ai , aj such that mij > 0 is connected by a straight-line

segment. In this way, a further partition of ā„¦ with an interesting property

results.

Theorem 6.4 If all Voronoi vertices are regular, then the dual graph coin-

cides with the set of edges of a triangulation of the convex hull of the given

point set.

This triangulation is called a Delaunay triangulation.

If among the Voronoi vertices there are degenerate ones, then a tri-

angulation can be obtained from the dual graph by a subsequent local

triangulation of the remaining m-polygons (m ā„ 4). A Delaunay triangula-

tion has the interesting property that two interior angles subtended by any

given edge sum to no more than Ļ. In this respect Delaunay triangulations

satisfy the ļ¬rst part of the angle condition formulated in Section 3.9 for

the maximum principle in ļ¬nite element methods.

Therefore, if ā„¦ is convex, then we automatically get a triangulation to-

gether with the Voronoi diagram. In the case of a nonconvex domain ā„¦,

certain modiļ¬cations could be required to achieve a correct triangulation.

.

264 6. Finite Volume Method

. .

.. This edge has to be

. removed from the

.. Delaunay triangulation.

Figure 6.7. Delaunay triangulation to the Voronoi diagram from Figure 6.5.

The implication

ā’

Voronoi diagram Delaunay triangulation ,

which we have just discussed, suggests that we ask about the converse

statement. We do not want to answer it completely at this point, but we

give the following suļ¬cient condition.

Theorem 6.5 If a conforming triangulation of ā„¦ (in the sense of ļ¬nite

element methods) consists of nonobtuse triangles exclusively, then it is

a Delaunay triangulation, and the corresponding Voronoi diagram can be

constructed by means of the perpendicular bisectors of the trianglesā™ edges.

We mention that the centre of the circumcircle of a nonobtuse triangle is

located within the closure of that triangle.

In the analysis of the ļ¬nite volume method, the following relation is

important.

Lemma 6.6 Given a nonobtuse triangle K with vertices aik , k ā {1, 2, 3},

then for the corresponding parts ā„¦ik ,K := ā„¦ik ā© K of the control volumes

ā„¦ik , we have

1 1

|K| ā¤ |ā„¦ik ,K | ā¤ |K| , k ā {1, 2, 3} .

4 2

The Donald diagram

In contrast to the Voronoi diagram, where the construction starts from a

given point set, the starting point here is a triangulation Th of ā„¦, which is

allowed to contain obtuse triangles.

Again, let K be a triangle with vertices aik , k ā {1, 2, 3}. We deļ¬ne

ā„¦ik ,K := x ā K Ī»j (x) < Ī»k (x), j = k ,

where Ī»k denote the barycentric coordinates with respect to aik (cf. (3.51)).

Obviously, the barycentre satisļ¬es aS = 1 (ai1 + ai2 + ai3 ), and (see, for

3

comparison, Lemma 6.6)

3 |ā„¦ik ,K | = |K| , k ā {1, 2, 3} . (6.4)

This relation is a simple consequence of the geometric interpretation of

the barycentric coordinates as area coordinates given in Section 3.3. The

6.2. Finite Volume Method on Triangular Grids 265

.

ai2

ā„¦i2 ,K

.

. aS

ā„¦i1 ,K ai1

ā„¦ i3 ,K

.

ai3

Figure 6.8. The subdomains ā„¦ik ,K .

required control volumes are deļ¬ned as follows (see Figure 6.8):

i ā Ī.

ā„¦i := int ā„¦i,K ,

K:ā‚K ai

The family {ā„¦i }iāĪ is called a Donald diagram.

The quantities Ī“ij , mij , and Īi are deļ¬ned similarly as in the case of

the Voronoi diagram. We mention that the boundary pieces Ī“ij are not

necessarily straight, but polygonal in general.

6.2.2 Finite Volume Discretization

The model under consideration is a special case of equation (6.1). Instead

of the matrix-valued diļ¬usion coeļ¬cient K we will take a scalar coeļ¬cient

k : ā„¦ ā’ R, that is, K = kI. Moreover, homogeneous Dirichlet boundary

conditions are to be satisļ¬ed. So the boundary value problem reads as

follows:

ā’ā Ā· (k āu ā’ c u) + r u =f in ā„¦ ,

(6.5)

u =0 on ā‚ā„¦ ,

with k, r, f : ā„¦ ā’ R, c : ā„¦ ā’ R2 .

The Case of the Voronoi Diagram

Let the domain ā„¦ be partioned by a Voronoi diagram and the correspond-

ing Delaunay triangulation. Due to the homogeneous Dirichlet boundary

conditions, it is suļ¬cient to consider only those control volumes ā„¦i that

are associated with inner nodes ai ā ā„¦. Therefore, we denote the set of

indices of all inner nodes by

Ī := i ā Ī ai ā ā„¦ .

In the ļ¬rst step, the diļ¬erential equation (6.5) is integrated over the single

control volumes ā„¦i :

ā’ ā Ā· (k āu ā’ c u) dx + i ā Ī.

r u dx = f dx , (6.6)

ā„¦i ā„¦i ā„¦i

266 6. Finite Volume Method

The application of Gaussā™s divergence theorem to the ļ¬rst integral of the

left-hand side of (6.6) yields

ā Ā· (k āu ā’ c u) dx = Ī½ Ā· (k āu ā’ c u) dĻ .

ā„¦i ā‚ā„¦i

Due to ā‚ā„¦i = āŖjāĪi Ī“ij (cf. Figure 6.9), it follows that

ā Ā· (k āu ā’ c u) dx = Ī½ij Ā· (k āu ā’ c u) dĻ ,

ā„¦i Ī“ij

jāĪi

where Ī½ij is the (constant) outer unit normal to Ī“ij (with respect to ā„¦i ).

In the next step we approximate the line integrals over Ī“ij .

.

.

. ai

Ī“ij

aj

.

Figure 6.9. The edge Ī“ij .

First, the coeļ¬cients k and Ī½ij Ā· c are approximated on Ī“ij by constants

Āµij > 0, respectively Ī³ij :

k|Ī“ij ā Āµij = const > 0 , Ī½ij Ā· c|Ī“ij ā Ī³ij = const .

In the simplest case, the approximation can be realized by the correspond-

ing value at the midpoint aĪ“ij of the straight-line segment Ī“ij . A better

choice is

ļ£±

ļ£²1 Ī½ij Ā· c dĻ , mij > 0 ,

mij Ī“ij

Ī³ij := (6.7)

ļ£³

Ī½ij Ā· c(aĪ“ij ) , mij = 0 .

We thus obtain

ā Ā· (k āu ā’ c u) dx ā [Āµij (Ī½ij Ā· āu) ā’ Ī³ij u] dĻ .

ā„¦i Ī“ij

jāĪi

The normal derivatives are approximated by diļ¬erence quotients; that is,

u(aj ) ā’ u(ai )

Ī½ij Ā· āu ā with dij := |ai ā’ aj | .

dij

This formula is exact for such functions that are linear along the straight-

line segment between the points ai , aj . So it remains to approximate the

integral of u over Ī“ij . For this, a convex combination of the values of u at

6.2. Finite Volume Method on Triangular Grids 267

the nodes ai and aj is taken:

u|Ī“ij ā rij u(ai ) + (1 ā’ rij ) u(aj ) ,

where rij ā [0, 1] is a parameter to be deļ¬ned subsequently. In general, rij

depends on Āµij , Ī³ij , and dij .

Collecting all the above approximations, we arrive at the following

relation:

ā Ā· (k āu ā’ c u) dx

ā„¦i

u(aj ) ā’ u(ai )

ā ā’ Ī³ij [rij u(ai ) + (1 ā’ rij ) u(aj )]

Āµij mij .

dij

jāĪi

To approximate the remaining integrals from (6.6), the following formulas

are used:

r u dx ā with mi := |ā„¦i | ,

r(ai ) u(ai ) mi =: ri u(ai ) mi ,

ā„¦i

ā

f dx f (ai ) mi =: fi mi .

ā„¦i

Instead of ri := r(ai ) or fi := f (ai ), the approximations

1 1

ri := r dx respectively fi := f dx (6.8)

mi mi

ā„¦i ā„¦i

can also be used. Denoting the unknown approximate values for u(ai ) by

ui , we obtain the following linear system of equations:

ui ā’ uj

+ Ī³ij [rij ui + (1 ā’ rij ) uj ] mij + ri ui mi

Āµij

dij (6.9)

ńņš. 10 |