<<

. 11
( 16)



>>

j∈Λi
= 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

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

<<

. 11
( 16)



>>