<< . .

. 83
( : 95)



. . >>

Galerkin method in the H1 -norm, while the accuracy of the UP scheme in the
L2 -norm deteriorates dramatically because of the e¬ect of the arti¬cial viscosity
which is O(h). Conversely, the SG converges quadratically since the introduced
numerical viscosity is in this case O(h2 ) as h tends to zero. •
572 12. Two-Point Boundary Value Problems

102


101


100
||u ’ u h || 1
H (0,1)
’1
10


10’2

||u ’ u h || 2
L (0,1)
10’3


10’4 ’3
10’2 10’1
10


FIGURE 12.11. Convergence analysis for an advection-di¬usion problem


12.6 A Quick Glance to the Two-Dimensional Case
The game that we want to play is to extend (in a few pages) the basic ideas
illustrated so far to the two-dimensional case. The obvious generalization of
problem (12.1)-(12.2) is the celebrated Poisson problem with homogeneous
Dirichlet boundary condition

’ u=f in „¦,
(12.90)
u=0 on ‚„¦,

where u = ‚ 2 u/‚x2 + ‚ 2 u/‚y 2 is the Laplace operator and „¦ is a two-
dimensional bounded domain whose boundary is ‚„¦. If we allow „¦ to be
the unit square „¦ = (0, 1)2 , the ¬nite di¬erence approximation of (12.90)
that mimics (12.10) is

Lh uh (xi,j ) = f (xi,j ) for i, j = 1, . . . , N ’ 1,
(12.91)
uh (xi,j ) = 0 if i = 0 or N, j = 0 or N,

where xi,j = (ih, jh) (h = 1/N > 0) are the grid points and uh is a grid
function. Finally, Lh denotes any consistent approximation of the operator
L = ’ . The classical choice is

1
(4ui,j ’ ui+1,j ’ ui’1,j ’ ui,j+1 ’ ui,j’1 ) ,
Lh uh (xi,j ) = (12.92)
h2

where ui,j = uh (xi,j ), which amounts to adopting the second-order centred
discretization of the second derivative (10.65) in both x and y directions
(see Figure 12.12, left).
The resulting scheme is known as the ¬ve-point discretization of the
Laplacian. It is an easy (and useful) exercise for the reader to check that
the associated matrix Afd has (N ’ 1)2 rows and is pentadiagonal, with the
12.6 A Quick Glance to the Two-Dimensional Case 573




xi+1,j
xij xi,j+1
xi,j’1

xi’1,j



FIGURE 12.12. Left: ¬nite di¬erence grid and stencil for a squared domain.
Right: ¬nite element triangulation of a region around a hole


i-th row given by
±
 4 if j = i,

1
’1 if j = i ’ N ’ 1, i ’ 1, i + 1, i + N + 1,
(afd )ij = 2 (12.93)
h 

0 otherwise.

Moreover, Afd is symmetric positive de¬nite and it is also an M-matrix (see
Exercise 14). As expected, the consistency error associated with (12.92) is
second-order with respect to h, and the same holds for the discretization
error u ’ uh h,∞ of the method. More general boundary conditions than
the one considered in (12.90) can be dealt with by properly extending to
the two-dimensional case the mirror imaging technique described in Section
12.2.3 and in Exercise 10 (for a thorough discussion of this subject, see, e.g.,
[Smi85]).

The extension of the Galerkin method is (formally speaking) even more
straightforward and actually is still readable as (12.43) with, however, the
implicit understanding that both the function space Vh and the bilinear
form a(·, ·) be adapted to the problem at hand. The ¬nite element method
corresponds to taking

Vh = vh ∈ C 0 („¦) : vh |T ∈ Pk (T ) ∀T ∈ Th , vh |‚„¦ = 0 , (12.94)

where Th denotes here a triangulation of the domain „¦ as previously intro-
duced in Section 8.5.2, while Pk (k ≥ 1) is the space of piecewise polyno-
mials de¬ned in (8.35). Note that „¦ needs not to be a rectangular domain
(see Figure 12.12, right).
574 12. Two-Point Boundary Value Problems

As of the bilinear form a(·, ·), the same kind of mathematical manipula-
tions performed in Section 12.4.1 lead to

∇uh · ∇vh dxdy,
a(uh , vh ) =
„¦

where we have used the following Green™s formula that generalizes the
formula of integration by parts

’ ∇u · ∇v dxdy ’ ∇u · n v dγ,
u v dxdy = (12.95)
„¦ „¦ ‚„¦

for any u, v smooth enough and where n is the outward normal unit vector
on ‚„¦ (see Exercise 15).
The error analysis for the two-dimensional ¬nite element approximation
of (12.90) can still be performed through the combined use of Ce`™s lemma
a
and interpolation error estimates as in Section 12.4.5 and is summarized in
the following result, which is the two-dimensional counterpart of Property
12.1 (for its proof we refer, e.g., to [QV94], Theorem 6.2.1).

Property 12.2 Let u ∈ H1 („¦) be the exact solution of (12.90) and uh ∈ Vh
0
be its ¬nite element approximation using continuous piecewise polynomials
of degree k ≥ 1. Assume also that u ∈ Hs („¦) for some s ≥ 2. Then the
following error estimate holds
M
u ’ uh ¤ Chl u (12.96)
H1 („¦) Hl+1 („¦)
0 ±0
where l = min(k, s ’ 1). Under the same assumptions, one can also prove
that

u ’ uh ¤ Chl+1 u Hl+1 („¦) . (12.97)
L2 („¦)

We notice that, for any integer s ≥ 0, the Sobolev space Hs („¦) introduced
above is de¬ned as the space of functions with the ¬rst s partial derivatives
(in the distributional sense) belonging to L2 („¦). Moreover, H1 („¦) is the
0
1
space of functions of H („¦) such that u = 0 on ‚„¦. The precise mathemat-
ical meaning of this latter statement has to be carefully addressed since,
for instance, a function belonging to H1 („¦) does not necessarily mean to
0
be continuous everywhere. For a comprehensive presentation and analysis
of Sobolev spaces we refer to [Ada75] and [LM68].
Following the same procedure as in Section 12.4.3, we can write the ¬nite
element solution uh as
N
uh (x, y) = uj •j (x, y)
j=1
12.7 Applications 575

N
where {•j }j=1 is a basis for Vh . An example of such a basis in the case
k = 1 is provided by the so-called hat functions introduced in Section 8.5.2
(see Figure 8.7, right). The Galerkin ¬nite element method leads to the
solution of the linear system Afe u = f , where (afe )ij = a(•j , •i ).
Exactly as happens in the one-dimensional case, the matrix Afe is sym-
metric positive de¬nite and, in general, sparse, the sparsity pattern being
strongly dependent on the topology of Th and the numbering of its nodes.
Moreover, the spectral condition number of Afe is still O(h’2 ), which im-
plies that solving iteratively the linear system demands for the use of a
suitable preconditioner (see Section 4.3.2). If, instead, a direct solver is
used, one should resort to a suitable renumbering procedure, as explained
in Section 3.9.



12.7 Applications
In this section we employ the ¬nite element method for the numerical
approximation of two problems arising in ¬‚uid mechanics and in the sim-
ulation of biological systems.


12.7.1 Lubrication of a Slider
Let us consider a rigid slider moving in the direction x along a physical
support from which it is separated by a thin layer of a viscous ¬‚uid (which
is the lubricant). Suppose that the slider, of length L, moves with a ve-
locity U with respect to the plane support which is supposed to have an
in¬nite length. The surface of the slider that is faced towards the support
is described by the function s (see Figure 12.13, left).
Denoting by µ the viscosity of the lubricant, the pressure p acting on the
slider can be modeled by the following Dirichlet problem
±
 s3

’ = ’(U s) x ∈ (0, L),
p


 p(0) = 0, p(L) = 0.

Assume in the numerical computations that we are working with a conver-
gent-divergent slider of unit length, whose surface is s(x) = 1’3/2x+9/8x2
with µ = 1.
Figure 12.13 (right) shows the solution obtained using linear and quadratic
¬nite elements with an uniform grid size h = 0.2. The linear system has
been solved by the nonpreconditioned conjugate gradient method. To re-
duce the Euclidean norm of the residual below 10’10 , 4 iterations are
needed in the case of linear ¬nite elements while 9 are required for quadratic
¬nite elements.
576 12. Two-Point Boundary Value Problems


L
0.09

U 0.08


0.07


0.06


0.05


0.04


s(x) 0.03


0.02



x 0.01


0


’0.01
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1




FIGURE 12.13. Left: geometrical parameters of the slider considered in Section
12.7.1; right: pressure on a converging-diverging slider. The solid line denotes
the solution obtained used linear ¬nite elements (the symbol —¦ denotes the nodal
values), while the dashed line denotes the solution obtained using quadratic ¬nite
elements

Table 12.2 reports a numerical study of the condition number K2 (Afe ) as
a function of h. In the case of linear ¬nite elements we have denoted the
matrix by A1 , while A2 is the corresponding matrix for quadratic elements.
Here we assume that the condition number approaches h’p as h tends to
zero; the numbers pi are the estimated values of p. As can be seen, in
both cases the condition number grows like h’2 , however, for every ¬xed
h, K2 (A2 ) is much bigger than K2 (A1 ).

h K2 (A1 ) p1 K2 (A2 ) p2
0.10000 63.951 “ 455.24
0.05000 348.21 2.444 2225.7 2.28
0.02500 1703.7 2.290 10173.3 2.19
0.01250 7744.6 2.184 44329.6 2.12
0.00625 33579 2.116 187195.2 2.07
TABLE 12.2. Condition number of the sti¬ness matrix for linear and quadratic
¬nite elements




12.7.2 Vertical Distribution of Spore Concentration over
Wide Regions
In this section we are concerned with di¬usion and transport processes of
spores in the air, such as the endspores of bacteria and the pollen of ¬‚ow-
ering plants. In particular, we study the vertical concentration distribution
of spores and pollen grains over a wide area. These spores, in addition to
settling under the in¬‚uence of gravity, di¬use passively in the atmosphere.
12.7 Applications 577

The basic model assumes the di¬usivity ν and the settling velocity β to
be given constants, the averaging procedure incorporating various physical
processes such as small-scale convection and horizontal advection-di¬usion
which can be neglected over a wide horizontal area. Denoting by x ≥ 0 the
vertical direction, the steady-state distribution u(x) of the spore concen-
tration is the solution of

’νu + βu = 0 0 < x < H,
(12.98)
’νu (H) + βu(H) = 0,
u(0) = u0 ,

where H is a ¬xed height at which we assume a vanishing Neumann condi-
tion for the advective-di¬usive ¬‚ux ’νu +βu (see Section 12.4.1). Realistic
values of the coe¬cients are ν = 10 m2 s’1 and β = ’0.03 ms’1 ; as for
u0 , a reference concentration of 1 pollen grain per m3 has been used in the
numerical experiments, while the height H has been set equal to 10 km.
The global P´clet number is therefore Pegl = 15.
e
The Galerkin ¬nite element method with piecewise linear ¬nite elements
has been used for the approximation of (12.98). Figure 12.14 (left) shows
the solution computed by running Program 94 on a uniform grid with step-
size h = H/10. The solution obtained using the (non stabilized) Galerkin
formulation (G) is denoted by the solid line, while the dash-dotted and
dashed lines refer to the Scharfetter-Gummel (SG) and upwind (UP) sta-
bilized methods. Spurious oscillations can be noticed in the G solution while
the one obtained using UP is clearly overdi¬used with respect to the SG
solution that is nodally exact. The local P´clet number is equal to 1.5 in
e
this case. Taking h = H/100 yields for the pure Galerkin scheme a stable
result as shown in Figure 12.14 (right) where the solutions furnished by G
(solid line) and UP (dashed line) are compared.

1
1.2
0.9

1 0.8

0.7
0.8

0.6
0.6
0.5

0.4 0.4

0.3
0.2
0.2
0
0.1

0
’0.2
0 200 400 600 800 1000 1200 1400 1600 1800
0 2000 4000 6000 8000 10000



FIGURE 12.14. Vertical concentration of spores: G, SG and UP solutions with
h = H/10 (left) and h = H/100 (right, where only the portion [0, 2000] is shown).
The x-axis represents the vertical coordinate
578 12. Two-Point Boundary Value Problems

12.8 Exercises
1. Consider the boundary value problem (12.1)-(12.2) with f (x) = 1/x. Using
(12.3) prove that u(x) = ’x log(x). This shows that u ∈ C 2 (0, 1) but u(0)
is not de¬ned and u , u do not exist at x = 0 (’: if f ∈ C 0 (0, 1), but not
f ∈ C 0 ([0, 1]), then u does not belong to C 0 ([0, 1])).
2. Prove that the matrix Afd introduced in (12.8) is an M-matrix.
[Hint: check that Afd x ≥ 0 ’ x ≥ 0. To do this, for any ± > 0 set Afd,± =
Afd +±In’1 . Then, compute w = Afd,± x and prove that min1¤i¤(n’1) wi ≥
0. Finally, since the matrix Afd,± is invertible, being symmetric and positive
de¬nite, and since the entries of A’1 are continuous functions of ± ≥ 0,
fd,±
’1
one concludes that Afd,± is a nonnegative matrix as ± ’ 0.]
0
3. Prove that (12.13) de¬nes a norm for Vh .
4. Prove (12.15) by induction on m.
5. Prove the estimate (12.23).
[Hint: for each internal node xj , j = 1, . . . , n ’ 1, integrate by parts (12.21)
to get

„h (xj )
® 
xj xj +h
1 
= ’u (xj ) ’ u (t)(xj ’ h ’ t)2 dt ’ u (t)(xj + h ’ t)2 dt» .
°
h2
xj ’h xj


Then, pass to the squares and sum „h (xj )2 for j = 1, . . . , n ’ 1. On noting
that (a + b + c)2 ¤ 3(a2 + b2 + c2 ), for any real numbers a, b, c, and applying
the Cauchy-Schwarz inequality yields the desired result.]
6. Prove that Gk (xj ) = hG(xj , xk ), where G is Green™s function introduced in
(12.4) and Gk is its corresponding discrete counterpart solution of (12.4).
[Solution: we prove the result by verifying that Lh G = hek . Indeed, for a
¬xed xk the function G(xk , s) is a straight line on the intervals [0, xk ] and
[xk , 1] so that Lh G = 0 at every node xl with l = 0, . . . , k ’ 1 and l =
k + 1, . . . , n + 1. Finally, a direct computation shows that (Lh G)(xk ) = 1/h
which concludes the proof.]
7. Let g = 1 and prove that Th g(xj ) = 1 xj (1 ’ xj ).
2
[Solution: use the de¬nition (12.25) with g(xk ) = 1, k = 1, . . . , n ’ 1 and
recall that Gk (xj ) = hG(xj , xk ) from the exercise above. Then
® 

<< . .

. 83
( : 95)



. . >>