= f i mi , i ∈ Λ .

This representation clearly indicates the a¬nity of the ¬nite volume method

to the ¬nite di¬erence method. However, for the subsequent analysis it is

more convenient to rewrite this system of equations in terms of a discrete

variational equality.

Multiplying the ith equation in (6.9) by arbitrary numbers vi ∈ R and

summing the results up over i ∈ Λ, we get

±

ui ’ uj

+ γij [rij ui + (1 ’ rij ) uj ] mij + ri ui mi

vi µij

dij

i∈Λ j∈Λi

= fi vi mi .

i∈Λ

Further, let Vh denote the space of continuous functions that are piecewise

linear over the (Delaunay) triangulation of „¦ and that vanish on ‚„¦. Then

the values ui and vi can be interpolated in Vh ; that is, there are unique

268 6. Finite Volume Method

uh , vh ∈ Vh such that uh (ai ) = ui , vh (ai ) = vi for all i ∈ Λ. The following

discrete bilinear forms on Vh — Vh can then be de¬ned:

mij

µij (ui ’ uj )

a0 (uh , vh ) := vi ,

h

dij

i∈Λ j∈Λi

[rij ui + (1 ’ rij ) uj ] γij mij ,

bh (uh , vh ) := vi

i∈Λ j∈Λi

dh (uh , vh ) := ri ui vi mi ,

i∈Λ

a0 (uh , vh )

ah (uh , vh ) := + bh (uh , vh ) + dh (uh , vh ) .

h

Finally, for two continuous functions v, w ∈ C(„¦), we set

w, v := wi vi mi ,

0,h

i∈Λ

where vi := v(ai ), wi := w(ai ).

Remark 6.7 ·, · 0,h is a scalar product on Vh . In particular, the following

norm can be introduced:

vh ∈ Vh .

vh := vh , vh 0,h , (6.10)

0,h

In (3.136) a discrete (L2 -) norm for a general ¬nite element space vh has

been de¬ned using the same notation. This multiple use seems to be accept-

able, since for regular triangulations both norms are equivalent uniformly

in h (see Remark 6.16 below).

Now the discrete variational formulation of the ¬nite volume method is

this:

Find uh ∈ Vh such that

for all vh ∈ Vh .

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

0,h

Up to now, the choice of the weighting parameters rij has remained open.

For this, two cases can be roughly distinguished:

(1) There exists a pair of indices (i, j) ∈ Λ — Λ such that µij |γij |dij .

|γij |dij .

(2) There is no such pair (i, j) with µij

In the second case, an appropriate choice is rij ≡ 1 . To some extent, this

2

can be seen as a generalization of the central di¬erence method to nonuni-

form grids. The ¬rst case corresponds to a locally convection-dominated

situation and requires a careful selection of the weighting parameters rij .

This will be explained in more detail in Section 9.3.

In general, the weighting parameters are of the following structure:

γij dij

rij = R , (6.12)

µij

6.2. Finite Volume Method on Triangular Grids 269

γ d

where R : R ’ [0, 1] is some function to be speci¬ed. The argument ij ijij

µ

is called the local P´clet number. Typical examples for this function R are

e

1

R(z) = [sign (z) + 1] , full upwinding,

2

(1 ’ „ )/2 , z < 0 , 2

„ (z) := max 0, 1 ’

R(z) = ,

(1 + „ )/2 , z ≥ 0 , |z|

1 z

R(z) = 1 ’ 1’ z , exponential upwinding .

e ’1

z

All these functions possess many common properties. For example, for all

z ∈ R,

[1 ’ R(z) ’ R(’z)] z

(P1) = 0,

R(z) ’ 1 z ≥

(P2) 0, (6.13)

2

1 ’ [1 ’ R(z)] z ≥

(P3) 0.

1

Note that the constant function R = satis¬es the conditions (P1) and

2

(P2) but not (P3).

The Case of the Donald Diagram

Let the domain „¦ be triangulated as in the ¬nite element method. Then,

following the explanations given in the second part of Section 6.2.1, the

corresponding Donald diagram can be created.

The discrete bilinear form in this case is de¬ned by

ah (uh , vh ) := k ∇uh , ∇vh + bh (uh , vh ) + dh (uh , vh ) ;

0

that is, the principal part of the di¬erential expression is discretized as in

the ¬nite element method, where bh , dh , and Vh are de¬ned as in the ¬rst

part of this section.

6.2.3 Comparison with the Finite Element Method

As we have already seen in Example 6.1, it may happen that a ¬nite volume

discretization coincides with a ¬nite di¬erence or ¬nite element discretiza-

tion. We also mention that the control volumes from that example are

exactly the Voronoi polygons to the grid points (i.e., to the nodes of the

triangulation).

Here we will consider this observation in more detail. By {•i }i∈Λ we de-

note the nodal basis of the space Vh of continuous, piecewise linear functions

on a conforming triangulation of the domain „¦.

Lemma 6.8 Let Th be a conforming triangulation of „¦ (in the sense of

¬nite element methods), all triangles of which are nonobtuse, and consider

the corresponding Voronoi diagram in accordance with Theorem 6.5. Then,

for an arbitrary triangle K ∈ Th with vertices ai , aj (i = j), the following

270 6. Finite Volume Method

relation holds:

mKij

∇•j · ∇•i dx = ’ ,

dij

K

where mK is the length of the segment of “ij that intersects K.

ij

Proof: Here we use some of the notation and the facts prepared at the

beginning of Section 3.9. In particular, ±K denotes the interior angle of K

ij

that is located in opposite the edge with vertices ai , aj . Next, the following

equality is an obvious fact from elementary geometry: 2 sin ±K mK = ij ij

K

cos ±ij dij . It remains to recall the relation

1

∇•j · ∇•i dx = ’ cot ±K

ij

2

K

2

from Lemma 3.47, and the statement immediately follows.

Corollary 6.9 Under the assumptions of Lemma 6.8, we have for k ≡ 1,

∇uh , ∇vh for all uh , vh ∈ Vh .

= a0 (uh , vh )

h

0

Proof: It is su¬cient to verify the relation for vh = •i and arbitrary

i ∈ Λ. First, we see that

∇uh , ∇•i ∇uh · ∇•i dx .

=

0

K

K‚supp•i

Furthermore,

∇uh · ∇•i dx ∇•j · ∇•i dx

= uj

K K

j:‚K aj

∇•i · ∇•i dx + ∇•j · ∇•i dx .

= ui uj

K K

j=i:‚K aj

Since

1= •j

j:‚K aj

over K, it follows that

∇•i = ’ ∇•j ; (6.14)

j=i:‚K aj

that is, by means of Lemma 6.8,

∇uh · ∇•i dx (uj ’ ui ) ∇•j · ∇•i dx

=

K K

j=i:‚K aj

6.2. Finite Volume Method on Triangular Grids 271

mKij

(ui ’ uj )

= . (6.15)

dij

j=i:‚K aj

Summing over all K ‚ supp •i , we get

mij

∇uh , ∇•i (ui ’ uj ) 2

= a0 (uh , •i ) .

= h

0

dij

j∈Λi

Remark 6.10 By a more sophisticated argumentation it can be shown

that the above corollary remains valid if the di¬usion coe¬cient k is con-

stant on all triangles K ∈ Th and if the approximation µij is chosen

according to

±

1 k|K mK + k|K mK

ij ij

k dσ = , mij > 0 ,

µij := (6.16)

mij “ij mij

0, mij = 0 ,

where K, K are both triangles sharing the vertices ai , aj .

Treatment of Matrix-valued Di¬usion Coe¬cients

Corollary 6.9 and Remark 6.10 are valid only in the spatial dimension d =

2. However, for more general control volumes, higher spatial dimensions,

or not necessarily scalar di¬usion coe¬cients, weaker statements can be

proven.

As an example, we will state the following fact. As a by-product, we also

obtain an idea for how to derive discretizations in the case of matrix-valued

di¬usion coe¬cients. For a better distinction between the elements K of

the triangulation and the di¬usion coe¬cient, we keep the notation k for

the di¬usion coe¬cient, even if k is allowed to be a matrix-valued function

temporarily.

Lemma 6.11 Let Th be a conforming triangulation of „¦, where in the

case of the Voronoi diagram it is additionally required that all triangles be

nonobtuse. Furthermore, assume that the di¬usion matrix k : „¦ ’ R2,2 is

constant on the single elements of Th . Then for any i ∈ Λ and K ∈ Th we

have

(k∇uh ) · ∇•i dx = ’ (k∇uh ) · ν dσ for all uh ∈ Vh ,

‚„¦i ©K

K

where {„¦i }i∈Λ is either a Voronoi or a Donald diagram and ν denotes the

outer unit normal with respect to „¦i .

Without di¬culties, the proof can be carried over from the proof of a related

result in [20, Lemma 6.1].

Now we will show how to use this fact to formulate discretizations for the

case of matrix-valued di¬usion coe¬cients. Namely, using relation (6.14),

272 6. Finite Volume Method

we easily see that

(k∇uh ) · ν dσ uj (k∇•j ) · ν dσ

=

‚„¦i ©K ‚„¦i ©K

j:‚K aj

(uj ’ ui ) (k∇•j ) · ν dσ .

=

‚„¦i ©K

j=i:‚K aj

Summing over all triangles that lie in the support of •i , we obtain by

Lemma 6.11 the relation

(k∇uh ) · ∇•i dx = (ui ’ uj ) (k∇•j ) · ν dσ . (6.17)

„¦ ‚„¦i

j∈Λi

With the de¬nition

±

dij (k∇•j ) · ν dσ , mij > 0 ,

µij := (6.18)

m

ij ‚„¦i

0, mij = 0 ,

it follows that

mij

(k∇uh ) · ∇•i dx = µij (ui ’ uj ) .

dij

„¦ j∈Λi

Note that, in the case of Voronoi diagrams, (6.16) is a special case of the

choice (6.18).

Consequently, in order to obtain a discretization for the case of a matrix-

valued di¬usion coe¬cient, it is su¬cient to replace in the bilinear form bh

and, if the Voronoi diagram is used, also in a0 , the terms involving µij

h

according to formula (6.18).

Implementation of the Finite Volume Method

In principle, the ¬nite volume method can be implemented in di¬erent

ways. If the linear system of equations is implemented in a node-orientated

manner (as in ¬nite di¬erence methods), the entries of the system matrix

Ah and the components of the right-hand side q h can be taken directly

from (6.9).

On the other hand, an element-orientated assembling is possible, too.

This approach is preferable, especially in the case where an existing ¬nite

element program will be extended by a ¬nite volume module. The idea of

how to do this is suggested by equation (6.17). Namely, for any triangle

K ∈ Th , the restricted bilinear form ah,K with the appropriate de¬nition

of µij according to (6.18) is de¬ned as follows:

ah,K (uh , vh ) :=

±

ui ’ uj

+ γij [rij ui + (1 ’ rij ) uj ] mij + ri ui mi ,

K K

vi µij

j=i:

dij

i∈Λ

‚K aj

6.2. Finite Volume Method on Triangular Grids 273

where mK := |„¦i © K|. Then the contribution of the triangle K to the

i

matrix entry (Ah )ij of the matrix Ah is equal to ah,K (•j , •i ). In the same

way, the right-hand side of (6.9) can be split elementwise.

6.2.4 Properties of the Discretization

Here we will give a short overview of basic properties of ¬nite volume

methods. For the sake of simplicity, we restrict ourselves to the case of a

constant scalar di¬usion coe¬cient k > 0. Then, in particular, it is useful

to set µij := k for all i ∈ Λ, j ∈ Λi .

Lemma 6.12 Suppose the approximations γij of νij · c|“ij satisfy γji =

’γij and the rij are de¬ned by (6.12) with a function R satisfying (P1).

Then we get for all uh , vh ∈ Vh ,

1

bh (uh , vh ) = ui vi γij mij

2

i∈Λ j∈Λi

1 1 1

rij ’ (ui ’ uj ) (vi ’ vj ) + (uj vi ’ ui vj ) γij mij .

+

2 2 2

i∈Λ j∈Λi

Proof: First, we observe that bh can be rewritten as follows:

1

vi (1 ’ rij ) uj ’ ’ rij ui γij mij

bh (uh , vh ) =

2

i∈Λ j∈Λi

(6.19)

1

+ ui vi γij mij .

2

i∈Λ j∈Λi

In the ¬rst term, we change the order of summation and rename the indices:

1

vj (1 ’ rji ) ui ’ ’ rji uj γji mji

bh (uh , vh ) =

2

i∈Λ j∈Λi

1

+ ui vi γij mij .

2

i∈Λ j∈Λi

Next we make use of the following relations, which easily result from dji =

dij and the assumptions on γij and rij :

1 1

(1 ’ rji ) γji = ’rij γij , ’ rji ’ rij

γji = γij .

2 2

So we get, due to mji = mij ,

1

vj ’rij ui ’ ’ rij uj γij mij

bh (uh , vh ) =

2

i∈Λ j∈Λi

1

+ ui vi γij mij .

2

i∈Λ j∈Λi

274 6. Finite Volume Method

Taking the arithmetic mean of both representations of bh , we arrive at

1

bh (uh , vh ) = ui vi γij mij

2

i∈Λ j∈Λi

1 1

(1 ’ rij ) uj vi ’ rij ui vj ’ ’ rji (ui vi + uj vj ) γij mij

+

2 2

i∈Λ j∈Λi

1 1

’ rij (uj vi + ui vj ’ ui vi ’ uj vj )

=

2 2

i∈Λ j∈Λi

1 1

(uj vi ’ ui vj ) γij mij +

+ ui vi γij mij .

2 2

i∈Λ j∈Λi

2

Corollary 6.13 Let c1 , c2 , ∇ · c ∈ C(„¦). Under the assumptions of

Lemma 6.12 and also assuming property (P2) for R, the bilinear form bh

satis¬es for all vh ∈ Vh the estimate

®

1

vi ° (γij ’ νij · c) dσ » . (6.20)

bh (vh , vh ) ≥ ∇ · c dx +

2

2 „¦i “ij

i∈Λ j∈Λi

Proof: Due to rij ’ 1 γij ≥ 0, because of property (P2) in (6.13), it

2

immediately follows that

1 1

bh (vh , vh ) ≥ 2 2

vi γij mij = vi γij mij .

2 2

i∈Λ j∈Λi i∈Λ j∈Λi

For the inner sum, we can write

γij mij = γij dσ

“ij

j∈Λi j∈Λi

νij · c dσ + (γij ’ νij · c) dσ .

=

“ij “ij

j∈Λi j∈Λi

The ¬rst term can be rewritten as an integral over the boundary of „¦i , i.e.,

νij · c dσ = ν · c dσ .

“ij ‚„¦i

j∈Λi

By Gauss™s divergence theorem, it follows that

ν · c dσ = ∇ · c dx .

‚„¦i „¦i

2

6.2. Finite Volume Method on Triangular Grids 275

Remark 6.14 If the approximations γij are chosen according to (6.7),

then γji = ’γij , and (6.20) simpli¬es to

1

bh (vh , vh ) ≥ ∇ · c dx .

2

vi

2 „¦i

i∈Λ

Using a similar argument as in the treatment of the term j∈Λi γij mij

in the proof of Corollary 6.13, the value dh (vh , vh ) can be represented as

follows:

2 2

dh (vh , vh ) = ri vi mi = vi ri dx

„¦i

i∈Λ i∈Λ

(ri ’ r) dx .

2 2

= vi r dx + vi (6.21)

„¦i „¦i

i∈Λ i∈Λ

The second term vanishes if the approximations ri are de¬ned as in (6.8).

Theorem 6.15 Let the rij be de¬ned by (6.12) with R satisfying (P1) and

(P2). Suppose k > 0, c1 , c2 , ∇ · c, r ∈ C(„¦), r + 1 ∇ · c ≥ r0 = const ≥ 0 on

2

„¦ and that the approximations γij , respectively ri , are chosen according to

(6.7), respectively (6.8). Under the assumptions of Lemma 6.8, we have for

all vh ∈ Vh ,

2

ah (vh , vh ) ≥ k ∇vh , ∇vh vi mi = k |vh |1 + r0 vh

2 2

+ r0 0,h ;

0

i∈Λ

that is, the bilinear form ah is Vh -elliptic uniformly with respect to h.

Proof: We start with the consideration of a0 (vh , vh ). Due to Corollary 6.9,

h

the relation

2

a0 (vh , vh ) = k ∇vh , ∇vh = k |vh |1

h 0

holds. Furthermore, by Remark 6.14 and equation (6.21), we have

1

bh (vh , vh ) + dh (vh , vh ) ≥ ∇·c+r dx ≥ r0

2 2

vi vi mi .

2

„¦i

i∈Λ i∈Λ

Since by de¬nition,

ah (vh , vh ) = a0 (vh , vh ) + bh (vh , vh ) + dh (vh , vh ) ,

h

2

both relations yield the assertion.

Remark 6.16 Let the family of triangulations (Th )h be regular. Then the

norms de¬ned in (3.136) and in (6.10) and also the norms · 0,h and

· 0 are equivalent on Vh uniformly with respect to h; i.e., there exist two

constants C1 , C2 > 0 independent of h such that

¤v ¤ C2 v for all v ∈ Vh .

C1 v 0 0,h 0

276 6. Finite Volume Method

Proof: Due to Theorem 3.43 (i) only the uniform equivalence of the dis-

crete L2 -norms has to be shown. Denoting such an equivalence by ∼, we

=

have for v ∈ Vh with vi := v(ai ) for i ∈ Λ,

« 1/2

1/2

¬ ·

=¬ |„¦i,K |·

|vi |2 mi |vi |2

K∈Th :

i∈Λ i∈Λ

K©„¦i =…

« 1/2

¬ ·

∼¬ |vi |2 |K|1/2 ·

=

K∈Th :

i∈Λ

K©„¦i =…

due to Lemma 6.6 or (6.4)

« 1/2

∼¬ ·

|K| |vi |2

=

K∈Th i:

ai ∈K

« 1/2

∼¬ ·

|vi |2

h2

= ,

K

K∈Th i:

ai ∈K

since due to the regularity of (Th )h there is a uniform lower bound for

the angles of K ∈ Th (see (3.93)) and thus a uniform upper bound on the

number of K ∈ Th such that K © „¦i = ….

2

Corollary 6.17 Under the assumptions of Theorem 6.15 and for a regular

family of triangulations (Th )h there exists a constant ± > 0 independent of

h such that

ah (vh , vh ) ≥ ± vh for all vh ∈ Vh .

2

1

Proof: By Remark 6.16 and Theorem 6.15,

2

ah (vh , vh ) ≥ k |vh |1 + r0 C1 vh

2 2

,

0

i.e., we can take ± := min{k; r0 C1 } . 2

2

Theorem 6.15 (or Corollary 6.17) asserts the stability of the method. It

is the fundamental result for the proof of an error estimate.

Theorem 6.18 Let {Th }h∈(0,h] be a regular family of conforming triangu-

¯

lations, where in the case of the Voronoi diagram it is additionally required

that all triangles be nonobtuse. Furthermore, suppose in (6.5) that k > 0,

c1 , c2 , ∇ · c, r ∈ C(„¦), r + 1 ∇ · c ≥ r0 = const > 0 on „¦, f ∈ C 1 („¦), and

2

that the approximations γij , respectively ri , are chosen according to (6.7),

6.2. Finite Volume Method on Triangular Grids 277

respectively (6.8). Let the rij be de¬ned by (6.12) with R satisfying (P1)

and (P2). If the exact solution u of (6.5) belongs to H 2 („¦) and uh ∈ Vh

denotes the solution of (6.11), then

u ’ uh ¤ C h[ u + |f |1,∞ ] ,

2

1

where the constant C > 0 is independent of h.

Proof: The proof rests on a similar idea to those in the proof and the

application of Strang™s ¬rst lemma (Theorem 3.38) in Section 3.6.

Denoting by Ih : C „¦ ’ Vh the interpolation operator de¬ned in (3.71)

¯

and setting vh := uh ’ Ih (u), we have

ah (vh , vh ) = ah (uh , vh ) ’ ah (Ih (u), vh )

= f, vh 0,h ’ ah (Ih (u), vh )

’ f dx ’ ah (Ih (u), vh ).

= f, vh vi f dx + vi

0,h

„¦i „¦i

i∈Λ i∈Λ

By the de¬nition of the discrete form f, vh 0,h and by the di¬erential

equation (6.5), considered as an equation in L2 („¦), we get

(fi ’ f ) dx + Lu dx ’ ah (Ih (u), vh ) ,

ah (vh , vh ) = vi vi

„¦i „¦i

i∈Λ i∈Λ

where Lu = ’∇ · (k ∇u ’ c u) + r u .

For f ∈ C 1 („¦) and the choice fi := f (ai ), it is easy to see that

|fi ’ f (x)| ¤ |f |1,∞ max hK ¤ C h|f |1,∞ for all x ∈ „¦i .

K:ai ∈K

So it follows that

(fi ’ f ) dx ¤ C h|f |1,∞ |vi |mi

vi

„¦i

i∈Λ i∈Λ

1/2 1/2

¤ C h|f |1,∞ 2

vi mi mi

i∈Λ i∈Λ

√

¤ |„¦|

¤ C h|f |1,∞ vh 0,h .

For the other choice of fi (see (6.8)), the same estimate is trivially satis¬ed.

The di¬cult part of the proof is to get an estimate of the consistency error

Lu dx ’ ah (Ih (u), vh ) .

vi

„¦i

i∈Λ

278 6. Finite Volume Method

This is very extensive, and so we will omit the details. A complete proof of

the following result is given in the paper [40]:

1/2

Lu dx ’ ah (Ih (u), vh ) ¤ C h u |vh |2 + vh 2

vi . (6.22)

2 1 0,h

„¦i

i∈Λ

Putting both estimates together and taking into consideration Remark 6.16,

we arrive at

1/2

ah (vh , vh ) ¤ C h [ u + |f |1,∞ ] |vh |2 + vh 2

2 1 0,h

¤ C h[ u + |f |1,∞ ] vh .

2 1

By Corollary 6.17, we conclude from this that

¤ C h[ u + |f |1,∞ ] .

vh 1 2

It remains to apply the triangle inequality and the standard interpolation

error estimate (cf. Theorem 3.29 with k = 1 or Theorem 3.35)

u ’ uh ¤ u ’ Ih (u) ¤ C h[ u + |f |1,∞ ] .

+ vh

1 1 1 2

2

We point out that the error measured in the H 1 -seminorm is of the same

order as for the ¬nite element method with linear ¬nite elements.

Now we will turn to the investigation of some interesting properties of

the method.

Global Conservativity

Here we consider the boundary value problem

’∇ · (k ∇u ’ c u) = f in „¦ ,

ν · (k ∇u ’ c u) = g on ‚„¦ .

Integrating the di¬erential equation over „¦, we conclude from Gauss™s

divergence theorem that

’ ∇ · (k ∇u ’ c u) dx = ’ ν · (k ∇u ’ c u) dσ = ’ g dσ ,

„¦ ‚„¦ ‚„¦

and hence

g dσ + f dx = 0 .

‚„¦ „¦

This is a necessary compatibility condition for the data describing the bal-

ance between the total ¬‚ow over the boundary and the distributed sources.

We will demonstrate that the discretization requires a discretized version of

this compatibility condition, which is called discrete global conservativity.

Therefore, we ¬rst have to de¬ne the discretization for the above type of

boundary conditions. Obviously, for inner control volumes „¦i (i ∈ Λ), there

6.2. Finite Volume Method on Triangular Grids 279

is no need for any modi¬cations. So we have to consider only the boundary

control volumes „¦i (i ∈ ‚Λ := Λ \ Λ).

In the case of the Voronoi diagram we have

’ ∇ · (k ∇u ’ c u) dx = ’ ν · (k ∇u ’ c u) dσ

„¦i ‚„¦i

’ ν · (k ∇u ’ c u) dσ ’ ν · (k ∇u ’ c u) dσ (6.23)

=

‚„¦i ©‚„¦

“ij

j∈Λi

’ ν · (k ∇u ’ c u) dσ ’

= g dσ .

‚„¦i ©‚„¦

“ij

j∈Λi

Since the line integrals over “ij can be approximated in the standard way,

we get the following equation:

ui ’ uj

+ γij [rij ui + (1 ’ rij ) uj ]

vi µij mij (6.24)

dij

j∈Λi

i∈Λ

’ vi g dσ = fi vi mi ,

‚„¦i ©‚„¦

i∈Λ i∈Λ

where the ansatz and test space Vh consists of all continuous functions over

¯

„¦ that are piecewise linear with respect to the underlying triangulation

(that is, in the boundary nodes no function values are prescribed).

It is again assumed that the rij are de¬ned by (6.12) with a function R

satisfying (P1) and γji = ’γij .

Obviously, the particular function ih :≡ 1 belongs to Vh . So we are

allowed to set vh = ih in the discretization. Then, repeating the above

symmetry argument (cf. the proof of Lemma 6.12), we get

mij mij

µij (ui ’ uj ) =’ µij (ui ’ uj ) ,

dij dij

i∈Λ j∈Λi i∈Λ j∈Λi

that is,

mij

µij (ui ’ uj ) = 0.

dij

i∈Λ j∈Λi

On the other hand, using the same argument, we have

[rij ui + (1 ’ rij ) uj ] γij mij

i∈Λ j∈Λi

[rji uj + (1 ’ rji ) ui ] γji mji

=

i∈Λ j∈Λi

=’ [(1 ’ rij ) uj + rij ui ] γij mij . (6.25)

i∈Λ j∈Λi

280 6. Finite Volume Method

Consequently, this term vanishes, too. Because of

vi g dσ = g dσ ,

‚„¦i ©‚„¦ ‚„¦

i∈Λ

it follows that

’ ≈

g dσ = fi vi mi = f i mi f dx . (6.26)

‚„¦ „¦

i∈Λ i∈Λ

This is the mentioned compatibility condition. It ensures the solvability of

the discrete system (6.24).

In the case of the Donald diagram, we obviously have

k∇uh , ∇vh = 0.

0

Since the proof of (6.25) does not depend on the particular type of the

control volumes, the property of discrete global conservativity in the sense

of (6.26) is satis¬ed for the Donald diagram, too.

Inverse Monotonicity

The so-called inverse monotonicity is a further important property of the

boundary value problem (6.5) that is inherited by the ¬nite volume dis-

cretization without any additional restrictive assumptions. Namely, it is

well known that under appropriate assumptions on the coe¬cients, the

solution u is nonnegative if the (continuous) right-hand side f in (6.5) is

nonnegative in „¦.

We will demonstrate that this remains true for the approximative solu-

tion uh . Only at this place is the property (P3) of the weighting function

R used; the preceding results are also valid for the simple case R(z) ≡ 1 .

2

There is a close relation to the maximum principles investigated in Sec-

tions 1.4 and 3.9. However, the result given here is weaker, and the proof

is based on a di¬erent technique.

Theorem 6.19 Let the assumptions of Theorem 6.15 be satis¬ed, but R

in (6.12) has to satisfy (P1)“(P3). Further, suppose that f ∈ C(„¦) and

f (x) ≥ 0 for all x ∈ „¦. Moreover, in the case of the Donald diagram, only

the weighting function R(z) = 1 [sign (z) + 1] is permitted.

2

Then

uh (x) ≥ 0 for all x ∈ „¦ .

Proof: We start with the case of the Voronoi diagram. Let uh be the

solution of (6.11) with f (x) ≥ 0 for all x ∈ „¦. Then we have the following

additive decomposition of uh :

uh = u+ ’ u’ , where u+ := max {0, uh } .

h h h

In general, u+ , u’ do not belong to Vh . So we interpolate them in Vh and

h h

set in (6.11) vh := Ih (u’ ), where Ih : C „¦ ’ Vh is the interpolation

¯

h

6.2. Finite Volume Method on Triangular Grids 281

operator (3.71). It follows that

= ah (uh , vh ) = ah Ih (u+ ), Ih (u’ ) ’ ah Ih (u’ ), Ih (u’ ) .

0 ¤ f, vh 0,h h h h h

By Theorem 6.15, we have

2

k Ih (u’ ) ¤ ah Ih (u’ ), Ih (u’ ) ¤ ah Ih (u+ ), Ih (u’ ) .

h h h h h

1

If we were able to show that ah Ih (u+ ), Ih (u’ ) ¤ 0, then the theorem

h h

would be proven, because this relation implies Ih (u’ ) 1 = 0, and from this

h

’

we immediately get uh = 0, and so uh = uh ≥ 0.

+

Since u+ u’ = 0 for all i ∈ Λ, it follows from (6.19) in the proof of

i i

Lemma 6.12 that

bh Ih (u+ ), Ih (u’ ) = (1 ’ rij ) u+ u’ γij mij . (6.27)

j i

h h

i∈Λ j∈Λi

Furthermore, obviously dh Ih (u+ ), Ih (u’ ) = 0 holds. Thus

h h

µij +

ah Ih (u+ ), Ih (u’ ) uj + γij (1 ’ rij ) u+ u’ mij

’

= j i

h h

dij

i∈Λ j∈Λi

µij γij dij

(1 ’ rij ) u+ u’ mij .

’ 1’

= j i

dij µij

i∈Λ j∈Λi

Due to 1 ’ [1 ’ R(z)] z ≥ 0 for all z ∈ R (cf. property (P3) in (6.13)) and

u+ u’ ≥ 0, it follows that

j i

ah Ih (u+ ), Ih (u’ ) ¤ 0 .

h h

So it remains to investigate the case of the Donald diagram. The function

R(z) = 1 [sign (z) + 1] has the property

2

1

[1 ’ R(z)] z = [1 ’ sign (z)] z ¤ 0 for all z ∈ R ,

2

that is (cf. (6.27)),

bh Ih (u+ ), Ih (u’ ) ¤ 0 .

h h

Taking u+ u’ = 0 into consideration, we get

i i

ah Ih (u+ ), Ih (u’ ) k ∇Ih (u+ ), ∇Ih (u’ )

¤

h h h h 0

u+ u’ ∇•j , ∇•i

=k .

j i 0

i∈Λ j∈Λi

Now Lemma 3.47 implies that

k

ah Ih (u+ ), Ih (u’ ) ¤ ’ u+ u’ cot ±K + cot ±K ,

ij ij

j i

h h

2

i∈Λ j∈Λi

where K and K are a pair of triangles sharing a common edge with vertices

ai , aj .

282 6. Finite Volume Method

Since all triangles are nonobtuse, we have cot ±K ≥ 0, cot ±K ≥ 0, and

ij ij

hence

ah Ih (u+ ), Ih (u’ ) ¤ 0 .

h h

2

Exercises

6.2 Suppose that the domain „¦ ‚ R2 can be triangulated by means of

equilateral triangles with edge length h > 0 in an admissible way.

(a) Give the shape of the control domains in the case of the Voronoi and

the Donald diagrams.

(b) Using the control domains from subproblem (a), discretize the Poisson

equation with homogeneous Dirichlet boundary conditions by means

of the ¬nite volume method.

1

6.3 Formulate an existence result for the weak solution in H0 („¦) of the

boundary value problem (6.5) similar to Theorem 3.12. In particular, what

form will condition (3.17) take?

6.4 Verify Remark 6.7; i.e., show that ·, · possesses the properties of

0,h

a scalar product on Vh .

6.5 Prove Remark 6.16 in detail.

6.6 Verify or disprove the properties (P1)“(P3) for the three weighting

functions given before (6.13) and for R ≡ 1 .

2

6.7 Let K be a nonobtuse triangle with the vertices a1 , a2 , a3 . The length

of the segments “K := “ij © K is denoted by mK , and dij is the length

ij ij

K

of the edge connecting ai with aj . Finally, ±ij is the interior angle of K

opposite that edge.

Demonstrate the following relation: 2mK = dij cot ±K .

ij ij

6.8

(a) Formulate problem (6.11) in terms of an algebraic system of type

(1.31).

(b) Show that for the resulting matrix Ah ∈ RM1 ,M1 , where M1 is the

number of elements of the index set Λ, the following relation is valid:

AT 1 ≥ 0 . Here, as in Section 1.4, 0, respectively 1, denotes a vector

h

of dimension M1 whose components are all equal to 0, respectively 1.

(This is nothing other than the property (1.32)(3)(i) except for the

transpose of Ah .)

7

Discretization Methods for Parabolic

Initial Boundary Value Problems

7.1 Problem Setting and Solution Concept

In this section initial boundary value problems for the linear case of the

di¬erential equation (0.33) are considered. We choose the form (3.12) to-

gether with the boundary conditions (3.18)“(3.20), which have already been

discussed in Section 0.5. In Section 3.2 conditions have been developed to

ensure a unique weak solution of the stationary boundary value problem. In

contrast to Chapter 3, the heterogeneities are now allowed also to depend

on time t, but for the sake of simplicity we do not do so for the coe¬cients

in the di¬erential equations and the boundary conditions, which covers

most of the applications, for example from Chapter 0. Also for the sake

of simplicity, we take the coe¬cient in front of the time derivative to be

constant and thus 1 by a proper scaling. From time to time we will restrict

attention to homogeneous Dirichlet boundary conditions for further ease of

exposition. Thus the problem reads as follows:

The domain „¦ is assumed to be a bounded Lipschitz domain and we

suppose that “1 , “2 , “3 form a disjoint decomposition of the boundary ‚„¦

(cf. (0.39)):

‚„¦ = “1 ∪ “2 ∪ “3 ,

where “3 is a closed subset of the boundary.

In the space-time cylinder QT = „¦ — (0, T ), T > 0, and its boundary

ST = ‚„¦ — (0, T ) there are given functions f : QT ’ R, g : ST ’ R,

g(x, t) = gi (x, t) for x ∈ “i , i = 1, 2, 3, and u0 : „¦ ’ R. The problem is to

284 7. Discretization of Parabolic Problems

¬nd a function u : QT ’ R such that

‚u

+ Lu = f in QT ,

‚t

(7.1)

Ru = g on ST ,

„¦ — {0} ,

u = u0 on

where Lv denotes the di¬erential expression for some function v : „¦ ’ R,

(Lv) (x) := ’∇ · (K(x) ∇v(x)) + c(x) · ∇v(x) + r(x)v(x) (7.2)

with su¬ciently smooth, time-independent coe¬cients

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

The boundary condition is expressed by the shorthand notation Ru = g,

which means, for a function ± : “2 ’ R on ‚„¦,

• Neumann boundary condition (cf. (0.41) or (0.36))

K∇u · ν = ‚νK u = g1 on “1 — (0, T ) , (7.3)

• mixed boundary condition (cf. (0.37))

K∇u · ν + ±u = ‚νK u + ±u = g2 on “2 — (0, T ) , (7.4)

• Dirichlet boundary condition (cf. (0.38))

on “3 — (0, T ) .

u = g3 (7.5)

Thus the stationary boundary problem considered so far reads

for x ∈ „¦ ,

Lu(x) = f (x) (7.6)

for x ∈ ‚„¦ .

Ru(x) = g(x)

It is to be expected that both for the analysis and the discretization there

are strong links between (7.6) and (7.1). The formulation (7.1) in particular

includes the heat equation (cf. (0.20))

‚u

’ ∇ · (K∇u) = f in QT , (7.7)

‚t

or for constant scalar coe¬cients in the form (cf. (0.19))

‚u

’ ∆u = f in QT (7.8)

‚t

with appropriate initial and boundary conditions.

Again as in Chapter 1, one of the simplest cases will be, for two space

dimensions (d = 2), the case of a rectangle „¦ = (0, a) — (0, b) or even the

case d = 1 (with „¦ = (0, a)), for which (7.8) further reduces to

‚2

‚u

’ 2 u = 0 in QT = (0, a) — (0, T ). (7.9)

‚t ‚x

For problem (7.1), the following typical analytical questions arise:

7.1. Problem Setting and Solution Concept 285

• existence of (classical) solutions,

• properties of the (classical) solutions,

• weaker concepts of the solution.

As in the case of elliptic boundary value problems, the theory of classical

solutions requires comparatively strong assumptions on the data of the

initial-boundary value problem. In particular, along the edge ‚„¦ — {0}

of the space-time cylinder initial and boundary conditions meet, so that

additional compatibility conditions have to be taken into account.

Representation of Solutions in a Special Case

To enhance the familiarity with the problem and for further comparison we

brie¬‚y sketch a method, named separation of variables, by which closed-

form solutions in the form of in¬nite series can be obtained for special cases.

Also in these cases, the representations are not meant to be a numerical

method (by its evaluation), but only serve as a theoretical tool.

We start with the case of homogeneous data, i.e., f = 0, gi = 0 (i =

1, 2, 3), so that the process is determined only by the initial data u0 .

We assume a solution of (7.1) to have the form u(x, t) = v(t)w(x) with

v(t) = 0, w(x) = 0. This leads to

’Lw(x)

v (t)

x ∈ „¦, t ∈ (0, T ) .

= , (7.10)

v(t) w(x)

Therefore, the expressions in (7.10) must be constant, for example, equal

to ’» for » ∈ R. Therefore,

v (t) = ’»v(t) , t ∈ (0, T ), (7.11)

which for the initial conditions v(0) = 1 has the solution

v(t) = e’»t .

Furthermore, w has to satisfy

Lw(x) = »w(x) , x ∈ „¦ ,

(7.12)

x ∈ ‚„¦ .

Rw(x) = 0 ,

Such a function w : „¦ ’ R, w = 0, is called an eigenfunction for the eigen-

¯

value » of the boundary value problem (7.6). If (wi , »i ), i = 1, . . . , N, are

eigenfunctions/values for (7.6), then because of the superposition principle,

the function

N

ci e’»i t wi (x)

u(x, t) := (7.13)

i=1

286 7. Discretization of Parabolic Problems

is a solution of the homogeneous initial-boundary value problem for the

initial value

N

u0 (x) := ci wi (x) , (7.14)

i=1

where the ci ∈ R are arbitrary. If there are in¬nitely many eigenfunc-

tions/values (wi , »i ) and if the sums in (7.13) and (7.14) converge in such

a way that also the in¬nite series possesses the derivatives appearing in

(7.6), then also

∞

ci e’»i t wi (x)

u(x, t) = (7.15)

i=1

is a solution to

∞

u0 (x) = ci wi (x) . (7.16)

i=1

For an inhomogeneous right-hand side of the form

N

f (x, t) = fi (t)wi (x) (7.17)

i=1

the solution representation can be extended to (variation of constants

formula)

t

N N

ci e’»i t wi (x) + fi (s)e’»i (t’s) ds wi (x) ,

u(x, t) := (7.18)

i=1 i=1 0

and at least formally the sum can be replaced by the in¬nite series. To

verify (7.18) it su¬ces to consider the case u0 = 0, for which we have

t

N N

fi (s)e’»i (t’s) ds »i wi (x)

fi (t)wi (x) ’

(‚t u)(x, t) =

i=1 i=1 0

t

N (7.19)

fi (s)e’»i (t’s) ds wi (x)

f (x, t) ’ L

=

i=1 0

f (x, t) ’ L(u)(x, t) .

=

From these solution representations we can conclude that initial data (and

thus also perturbances contained in it) and also the in¬‚uence of the right-

hand side act only exponentially damped if all eigenvalues are positive.

For d = 1, „¦ = (0, a) and Dirichlet boundary conditions we have the

eigenfunctions

π

ν ∈ N,

wν (x) = sin ν x , (7.20)

a

7.1. Problem Setting and Solution Concept 287

for the eigenvalues

νπ 2

»ν = . (7.21)

a

If the initial data u0 has the representation

∞

π

u0 (x) = ci sin ν x , (7.22)

a

ν=1

then for example for f = 0 the (formal) solution reads

∞

π

ci e’»ν t sin ν

u(x, t) = x . (7.23)

a

ν=1

The eigenfunctions wν are orthogonal with respect to the scalar product

·, · 0 in L2 („¦), since they satisfy

0 for ν = µ,

π π

· ·

sin ν , sin µ = (7.24)

a

a a for ν = µ,

0

2

which can by checked by means of well-known identities for the trigono-

metric functions.

Therefore (see below (7.57)),

u0 , w ν 0

ci = , (7.25)

wν , wν 0

which is called the Fourier coe¬cient in the Fourier expansion of u0 .

Of course, the (wν , »ν ) depend on the boundary conditions. For Neumann

boundary conditions in x = 0 and x = a we have

wν (x) = cos ν π x , ν = 0, 1, . . . ,

a

(7.26)

2

νπ

ν

»= , ν = 0, 1, . . . .

a

The occurrence of w0 = 1, »0 = 0 re¬‚ects the nontrivial solvability of the

pure Neumann problem (which therefore is excluded by the conditions of

Theorem 3.15).

For Lu = ’∆u and „¦ = (0, a) — (0, b), eigenfunctions and eigenvalues

can be derived from the one-dimensional case because of

˜

’∆(v ν (x)˜µ (y)) = ’v ν (x)˜µ (y) ’ v ν (x)˜µ (y) = (»ν + »µ )v ν (x)˜µ (y).

v v v v

Therefore, for „¦ = (0, a)—(0, b) one has to choose the eigenfunctions/values

(v ν , »ν ) (in x, on (0, a)) for the required boundary conditions at x = 0 and

v˜

x = a, and (˜µ , »µ ) (in y, on (0, b)) for the required boundary conditions

at y = 0, y = b.

For Dirichlet boundary conditions everywhere this leads to

π π

wνµ (x, y) = sin ν x sin µ x (7.27)

a b

288 7. Discretization of Parabolic Problems

for the eigenvalues

2 2

νπ µπ

»νµ = +

a b

2 2

and »νµ ’ ∞ for ν ’ ∞ or

(i.e., the smallest eigenvalue is π + π

a b

µ ’ ∞).

As a further concluding example we note the case

x = 0 or x = a : u(x, y) = 0 for y ∈ [0, b] ,

y = 0 : ∇u · ν(x, y) = ’‚2 u(x, y) = 0 for x ∈ (0, a) ,

y = b : ∇u · ν(x, y) = ‚2 u(x, y) = 0 for x ∈ (0, a) .

Eigenfunctions:

π π

wνµ (x, y) = sin ν x cos µ y , (7.28)

a b

ν = 1, 2, . . . , µ = 0, 1, 2, . . . .

Eigenvalues:

2 2

π π

νµ

» =ν +µ .

a b

A Sketch of the Theory of Weak Solutions

As in the study of the elliptic boundary value problems (3.12), (3.18)“

(3.20), for equation (7.1) a weak formulation can be given that reduces the

requirements with respect to the di¬erentiability properties of the solution.

The idea is to treat time and space variables in a di¬erent way:

• For ¬xed t ∈ (0, T ), the function x ’ u(x, t) is interpreted

(1)

as a parameter-dependent element u(t) of some space V whose

elements are functions of x ∈ „¦. An obvious choice is (see

Subsection 3.2.1, (I)) the space

V = {v ∈ H 1 („¦) : v = 0 on “3 }.

• In a next step, that is, for varying t, a function t ’ u(t) results

with values in the (function) space V.

(2) In addition to V, a further space H = L2 („¦) occurs, from which the

initial value u0 is taken and which contains V as a dense subspace.

A subspace V is called dense in H if the closure of V with respect to

the norm on H coincides with H.

(3) The time derivative is understood in a generalized sense; see (7.29).

(4) The generalized solution t ’ u(t) is sought as an element of a function

space, the elements of which are “function-valued” (cf. (1)).

De¬nition 7.1 Let X denote one of the spaces H or V (in particular, this

means that the elements of X are functions on „¦ ‚ Rd ).

7.1. Problem Setting and Solution Concept 289

(i) The space C l ([0, T ], X), l ∈ N0 , consists of all continuous functions

v : [0, T ] ’ X that have continuous derivatives up to the order l on

[0, T ] with the norm

l

v (i) (t)

sup .

X

i=0 t∈(0,T )

For the sake of simplicity, the notation C([0, T ], X) := C 0 ([0, T ], X)

is used.

(ii) The space Lp ((0, T ), X) with 1 ¤ p ¤ ∞ consists of all functions on

(0, T ) — „¦ with the following properties:

v(t, ·) ∈ X for any t ∈ (0, T ), F ∈ Lp (0, T ) with F (t) := v(t, ·) .

X

Furthermore,

v := F .

Lp ((0,T ),X) Lp (0,T )

Remark 7.2 f ∈ L2 (QT ) ’ f ∈ L2 ((0, T ), H) .

Proof:

2

Basically, the proof is a consequence of Fubini™s theorem (see [1]).

Concerning the interpretation of the time derivative and of the weak

formulation, a comprehensive treatment is possible only within the frame-

work of the theory of distributions; thus a detailed explanation is beyond

the scope of this book. A short but mathematically rigorous introduction

can be found in the book [39, Chapter 23].

The basic idea consists in the following de¬nition:

A function u ∈ L2 ((0, T ), V ) is said to have a weak derivative w if the

following holds:

T T

∞

u(t) Ψ (t) dt = ’ for all Ψ ∈ C0 (0, T ) .

w(t) Ψ(t) dt (7.29)

0 0

du

Usually, this derivative w is denoted by or u .

dt

Remark 7.3 The integrals occurring above are to be understood as so-

called Bochner integrals and are extensions of the Lebesgue integral to

function-valued mappings. Therefore, equation (7.29) is an equality of

functions.

Before we give a weak formulation of (7.1), the following notion is worth

recalling:

u v dx for u, v ∈ H ,

u, v := (7.30)

0

„¦

[K∇u · ∇v + (c · ∇u + ru) v] dx+ ±uv dσ, u, v ∈ V. (7.31)

a(u, v) :=

„¦ “2

290 7. Discretization of Parabolic Problems

Let u0 ∈ H, f ∈ L2 ((0, T ), H), and in case of Dirichlet conditions we

restrict ourselves to the homogeneous case.

An element u ∈ L2 ((0, T ), V ) is called a weak solution of (7.1) if it has

a weak derivative du = u ∈ L2 ((0, T ), H) and the following holds

dt

d

u(t), v + a (u(t), v) = f (t), v + g1 (·, t)v dσ

0

dt “1

0

+ g2 (·, t)v dσ (7.32)

“2

for all v ∈ V and t ∈ (0, T ) ,

u(0) = u0 .

Due to u ∈ L2 ((0, T ), V ) and u ∈ L2 ((0, T ), H) , we also have u ∈

C ([0, T ], H) (see [12, p. 287]), so that the initial condition is meaningful in

the classical sense.

In what follows, the bilinear form a is assumed to be continuous on V —V

(see (3.2)) and V -elliptic (see (3.3)). The latter means that there exists a

number ± > 0 such that

a(v, v) ≥ ± v for all v ∈ V .

2

V

Lemma 7.4 Let a be a V -elliptic, continuous bilinear form, u0 ∈ H, and

f ∈ C ([0, T ], H), and suppose the considered boundary conditions are ho-

mogeneous. Then, for the solution u(t) of (7.32) the following estimate

holds:

t

’±t

e’±(t’s) ds

¤ u0 for all t ∈ (0, T ) .

u(t) e + f (s)

0 0 0

0

Proof: The following equations are valid almost everywhere in (0, T ).

Setting v = u(t), (7.32) reads as

u (t), u(t) + a(u(t), u(t)) = f (t), u(t) .

0 0

Using the relation

1d 1d d

2

u (t), u(t) = u(t), u(t) = u(t) = u(t) u(t)

0 0

0

0 0

2 dt 2 dt dt

and the V -ellipticity, it follows that

d

¤ f (t), u(t)

2

u(t) u(t) + ± u(t) .

0 0 V 0

dt

Now the simple inequality

¤ u(t)

u(t) 0 V

and the Cauchy“Schwarz inequality

¤ f (t)

f (t), u(t) u(t)

0 0

0

7.1. Problem Setting and Solution Concept 291

yield, after division by u(t) 0 , the estimate

d

¤ f (t)

u(t) + ± u(t) 0.

0 0

dt

Multiplying this relation by e±t , the relation

d ±t d

(e u(t) 0 ) = e±t + ±e±t u(t)

u(t) 0 0

dt dt

leads to

d ±t

(e u(t) 0 ) ¤ e±t f (t) .

0

dt

The integration over (0, t) results in

t

’ u(0) ¤

±t

e±s f (s)

e u(t) ds

0 0 0

0

for all t ∈ (0, T ). Multiplying this by e’±t and taking into consideration

the initial condition, we get the asserted relation

t

’±t ’±(t’s)

¤ u0

u(t) 0e + f (s) 0e ds .

0

0

2

As a consequence of this lemma, the uniqueness of the solution of (7.32)

is obtained.

Corollary 7.5 Let a be a V -elliptic, continuous bilinear form. Then there

exists at most one solution of (7.32).

Proof: Suppose there are two di¬erent solutions u1 (t), u2 (t) ∈ V. Then

the di¬erence v(t) := u1 (t) ’ u2 (t) solves a homogeneous problem of the

type (7.32) (i.e., with f = 0, u0 = 0). Lemma 7.4 immediately implies

v(t) 0 = 0 in [0, T ); that is, u1 (t) = u2 (t) for all t ∈ [0, T ). 2

There is a close relation between Lemma 7.4 and solution representations

such as (7.18) (with the sum being in¬nite). The eigenvalue problem (7.12)

is de¬ned as follows in its variational form (see also the end of Section 2.2):

De¬nition 7.6 A number » ∈ R is called an eigenvalue for the eigenvector

w ∈ V, w = 0, if

for all v ∈ V .

a(w, v) = » w, v 0