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

6µ

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

®