10 10

’6

10

’3

10

’8

10

’4

10

’10

10

’5

10

’12

10

’6

10 ’14

10

’1 0 1

10 10 10 10 20 30 40 50 60

FIGURE 12.6. Left: error curves for linear and quadratic elements. The dashed

and solid lines denote the H1 (0, L) and L2 (0, L) norms of the error in the case

k = 1, while the dot-line and dotted line denote the corresponding norms in the

case k = 2. Right: error curves for the spectral collocation method. The dashed

and solid lines denote the H1 (0, L) and L2 (0, L) norms of the error, respectively

12.4.7 Spectral Methods

It turns out that the spectral collocation method of Section 12.3 can be

regarded as a Galerkin method where the subspace is P0 and the integrals

n

are approximated by the Gauss-Lobatto quadrature formula. As a matter

of fact, the approximation of problem (12.38) is

¬nd un ∈ P0 : an (un , vn ) = (f, vn )n ∀vn ∈ P0 , (12.68)

n n

where an is the bilinear form that is obtained from the bilinear form a

by replacing exact integrals by the Gauss-Lobatto formula (12.37). For

problem (12.41) the associated bilinear form a was introduced in (12.44).

We would therefore obtain

an (un , vn ) = (±un , vn )n + (βun , vn )n + (γun , vn )n . (12.69)

This is no longer a Galerkin method, but is called a generalized Galerkin

approximation. Its analysis requires more care than for Galerkin methods,

560 12. Two-Point Boundary Value Problems

as already seen in Section 12.3 and depends on the Strang lemma (see

[QV94]). However, the same kind of error estimate (12.39) can be proven

in this case as well.

A further generalization combining the ¬nite element approach with

piecewise polynomials of high degree and Gauss-Lobatto integration on

each element yields the so-called spectral element method and the h ’ p

version of the ¬nite element method (here p stands for the polynomial de-

gree that we have denoted with n). In these cases convergence is achieved

letting simultaneously (or independently) h go to zero and p go to in¬nity.

(See, e.g., [BM92], [SS98]).

Example 12.2 We consider again the two-point boundary value problem (12.67)

and employ the spectral collocation method for its numerical approximation.

We show in Figure 12.6 (right) the error curves in the L2 (0, L) (solid line) and

H1 (0, L) (dashed line) norms as functions of the spectral degree n, with n = 4’k ,

k = 1, . . . , 5. Notice the high accuracy that is achieved, even when a small value

of n is used, due to the smoothness of the exact solution. Notice also that for

n ≥ 32 the accuracy is actually bounded by the e¬ect of rounding errors. •

12.5 Advection-Di¬usion Equations

Boundary value problems of the form (12.41) are used to describe processes

of di¬usion, advection and absorption (or reaction) of a certain quantity

which is identi¬ed with u(x). The term ’(±u ) is responsible for the di¬u-

sion, βu for the advection (or transport), γu for the absorption (if γ > 0).

In this section we focus on the case where ± is small compared with β (or

γ). In these cases, the Galerkin method that we introduced earlier might

be unsuitable for providing accurate numerical results. A heuristic explana-

tion can be drawn from the inequality (12.56), noticing that in this case the

constant M/±0 can be very large, hence the error estimate can be mean-

ingless unless h is much smaller than (M/±0 )’1 . For instance, if ± = µ,

γ = 0 and β = const 1, then ±0 = µ and M = µ + CP β. Similarly, if

2

± = µ, β = 0 and γ = const 1 then ±0 = µ and M = µ + CP γ.

To keep our analysis at the simplest possible level, we will consider the

following elementary two-point boundary value problem

’µu + βu = 0, 0 < x < 1,

(12.70)

u(0) = 0, u(1) = 1,

where µ and β are two positive constants such that µ/β 1. Despite

its simplicity, (12.70) provides an interesting paradigm of an advection-

di¬usion problem in which advection dominates di¬usion.

12.5 Advection-Di¬usion Equations 561

We de¬ne the global P´clet number as

e

|β|L

Pegl = (12.71)

2µ

where L is the size of the domain (equal to 1 in our case). The global P´clet

e

number measures the dominance of the advective term over the di¬usive

one.

Let us ¬rst compute the exact solution of problem (12.70). The charac-

teristic equation associated to the di¬erential equation is ’µ»2 + β» = 0

and admits the roots »1 = 0 and »2 = β/µ. Then

β

u(x) = C1 e»1 x + C2 e»2 x = C1 + C2 e µ x ,

where C1 and C2 are arbitrary constants. Imposing the boundary conditions

yields C1 = ’1/(eβ/µ ’ 1) = ’C2 , therefore

u(x) = (exp(βx/µ) ’ 1) / (exp(β/µ) ’ 1) .

If β/µ 1 we can expand the exponentials up to ¬rst order obtaining

β β

x + · · · ’ 1)/(1 + + · · · ’ 1)

u(x) = (1 + (β x/µ)/(β/µ) = x,

µ µ

thus the solution is close to the solution of the limit problem ’µu = 0,

which is a straight line interpolating the boundary data.

However, if β/µ 1 the exponentials attain big values so that

exp(β/µx) β

= exp ’ (1 ’ x) .

u(x)

exp(β/µ) µ

Since the exponent is big and negative the solution is almost equal to zero

everywhere unless a small neighborhood of the point x = 1 where the

term 1 ’ x becomes very small and the solution joins the value 1 with an

exponential behaviour. The width of the neighbourhood is of the order of

µ/β and thus it is quite small: in such an event, we say that the solution

exhibits a boundary layer of width O (µ/β) at x = 1.

12.5.1 Galerkin Finite Element Approximation

Let us discretize problem (12.70) using the Galerkin ¬nite element method

introduced in Section 12.4.5 with k = 1 (piecewise linear ¬nite elements).

The approximation to the problem is: ¬nd uh ∈ Xh such that

1

±

a(uh , vh ) = 0 1,0

∀vh ∈ Xh ,

(12.72)

u (0) = 0, u (1) = 1,

h h

562 12. Two-Point Boundary Value Problems

1,0

1

where the ¬nite element spaces Xh and Xh have been introduced in (8.22)

and (12.57) and the bilinear form a(·, ·) is

1

a(uh , vh ) = (µuh vh + βuh vh ) dx. (12.73)

0

Remark 12.5 (Advection-di¬usion problems in conservation form)

Sometimes, the advection-di¬usion problem (12.70) is written in the follow-

ing conservation form

’(J(u)) = 0, 0 < x < 1,

(12.74)

u(0) = 0, u(1) = 1,

where J(u) = µu ’βu is the ¬‚ux (already introduced in the ¬nite di¬erence

context in Section 12.2.3), µ and β are given functions with µ(x) ≥ µ0 > 0

for all x ∈ [0, 1]. The Galerkin approximation of (12.74) using piecewise

linear ¬nite elements reads: ¬nd uh ∈ Xh such that

1

1,0

∀vh ∈ Xh

b(uh , vh ) = 0,

1

(µuh ’ βuh )vh dx. The bilinear form b(·, ·) coincides

where b(uh , vh ) =

0

with the corresponding one in (12.73) when µ and β are constant.

Taking vh as a test function the generic basis function •i , (12.72) yields

1 1

i = 1, . . . , n ’ 1.

µuh •i dx + βuh •i dx = 0

0 0

n

Setting uh (x) = j=0 uj •j (x), and noting that supp(•i ) = [xi’1 , xi+1 ] the

above integral, for i = 1, . . . , n ’ 1, reduces to

®

xi+1 xi+1

xi

µ °ui’1 •i •i+1 dx»

2

•i’1 •i dx + ui (•i ) dx + ui+1

®

xi’1 xi’1 xi

xi+1 xi+1

xi

+β °ui’1 •i+1 •i dx» = 0.

•i’1 •i dx + ui •i •i dx + ui+1

xi’1 xi’1 xi

Assuming a uniform partition of [0, 1] with xi = xi’1 + h for i = 1, . . . , n,

h = 1/n, and noting that •j (x) = h if xj’1 ¤ x ¤ xj , •j (x) = ’ h if

1 1

xj ¤ x ¤ xj+1 , we deduce that

µ 1

(’ui’1 + 2ui ’ ui+1 ) + β (ui+1 ’ ui’1 ) = 0, i = 1, . . . , n ’ 1.

h 2

(12.75)

12.5 Advection-Di¬usion Equations 563

Multiplying by h/µ and de¬ning the local P´clet number to be

e

|β|h

Pe =

2µ

we ¬nally obtain

(Pe ’ 1) ui+1 + 2ui ’ (Pe + 1) ui’1 = 0. i = 1, . . . , n ’ 1. (12.76)

This is a linear di¬erence equation which admits a solution of the form

ui = A1 ρi + A2 ρi for suitable constants A1 , A2 (see Section 11.4), where

1 2

ρ1 and ρ2 are the two roots of the following characteristic equation

(Pe ’ 1) ρ2 + 2ρ ’ (Pe + 1) = 0.

±

Thus

√ 1 + Pe

’1 ± 1 + Pe2 ’ 1 ,

1 ’ Pe

ρ1,2 = =

1.

Pe ’ 1

Imposing the boundary conditions at x = 0 and x = 1 gives

n

A1 = 1/(1 ’ ), A2 = ’A1 ,

1+Pe

1’Pe

so that the solution of (12.76) is

i n

1 + Pe 1 + Pe

1’ / 1’

ui = i = 0, . . . , n.

1 ’ Pe 1 ’ Pe

We notice that if Pe > 1, a power with a negative base appears at the

numerator which gives rise to an oscillating solution. This is clearly visible

in Figure 12.7 where for several values of the local P´clet number, the

e

solution of (12.76) is compared with the exact solution (sampled at the

mesh nodes) corresponding to a value of the global P´clet equal to 50.

e

The simplest remedy for preventing the oscillations consists of choosing

a su¬ciently small mesh stepsize h in such a way that Pe < 1. However this

approach is often impractical: for example, if β = 1 and µ = 5 · 10’5 one

should take h < 10’4 which amounts to dividing [0, 1] into 10000 subin-

tervals, a strategy that becomes unfeasible when dealing with multidimen-

sional problems. Other strategies can be pursued, as will be addressed in

the next sections.

12.5.2 The Relationship between Finite Elements and Finite

Di¬erences; the Numerical Viscosity

To examine the behaviour of the ¬nite di¬erence method (FD) when applied

to the solution of advection-di¬usion problems and its relationship with the

564 12. Two-Point Boundary Value Problems

1

0.5

0

’0.5

0.75 0.8 0.85 0.9 0.95 1

FIGURE 12.7. Finite di¬erence solution of the advection-di¬usion problem

(12.70) (with Pegl = 50) for several values of the local P´clet number. Solid

e

line: exact solution, dot-dashed line: Pe = 2.63, dotted line: Pe = 1.28, dashed

line: Pe = 0.63

¬nite element method (FE), we again consider the one-dimensional problem

(12.70) with a uniform meshsize h.

To ensure that the local discretization error is of second order we approx-

imate u (xi ) and u (xi ), i = 1, . . . , n ’ 1, by the centred ¬nite di¬erences

(10.61) and (10.65) respectively (see Section 10.10.1). We obtain the fol-

lowing FD problem

±

’µ ui+1 ’ 2ui + ui’1 + β ui+1 ’ ui’1 = 0, i = 1, . . . , n ’ 1

h2 (12.77)

2h

u = 0, u = 1.

0 n

If we multiply by h for every i = 1, . . . , n ’ 1, we obtain exactly the same

equation (12.75) that was obtained using piecewise linear ¬nite elements.

The equivalence between FD and FE can be pro¬tably employed to de-

vise a cure for the oscillations arising in the approximate solution of (12.75)

when the local P´clet number is larger than 1. The important observation

e

here is that the instability in the FD solution is due to the fact that the

discretization scheme is a centred one. A possible remedy consists of ap-

proximating the ¬rst derivative by a one-sided ¬nite di¬erence according

to the direction of the transport ¬eld. Precisely, we use the backward dif-

ference if the convective coe¬cient β is positive and the forward di¬erence

otherwise. The resulting scheme when β > 0 is

ui+1 ’ 2ui + ui’1 ui ’ ui’1

’µ i = 1, . . . , n ’ 1

+β =0 (12.78)

h2 h

which, for µ = 0 reduces to ui = ui’1 and therefore yields the desired

constant solution of the limit problem βu = 0. This one-side discretization

12.5 Advection-Di¬usion Equations 565

of the ¬rst derivative is called upwind di¬erencing: the price to be paid for

the enhanced stability is a loss of accuracy since the upwind ¬nite di¬erence

introduces a local discretization error of O(h) and not of O(h2 ) as happens

using centred ¬nite di¬erences.

Noting that

ui ’ ui’1 ui+1 ’ ui’1 h ui+1 ’ 2ui + ui’1

’

= ,

h2

h 2h 2

the upwind ¬nite di¬erence can be interpreted as the sum of a centred ¬nite

di¬erence approximating the ¬rst derivative and of a term proportional to

the discretization of the second-order derivative. Consequently, (12.78) is

equivalent to

ui+1 ’ 2ui + ui’1 ui+1 ’ ui’1

’µh i = 1, . . . , n ’ 1 (12.79)

+β =0

h2 2h

where µh = µ(1 + Pe). This amounts to having replaced the di¬erential

equation (12.72) with the perturbed one

’µh u + βu = 0, (12.80)

then using centred ¬nite di¬erences to approximate both u and u . The

perturbation

βh

’µ Pe u = ’ u (12.81)

2

is called the numerical viscosity (or arti¬cial di¬usion). In Figure 12.8

a comparison between centred and upwinded discretizations of problem

(12.72) is shown.

More generally, we can resort to a centred scheme of the form (12.80)

with the following viscosity

µh = µ(1 + φ(Pe)) (12.82)

where φ is a suitable function of the local P´clet number satisfying

e

lim φ(t) = 0.

t’0+

Notice that when φ(t) = 0 for all t, one recovers the centred ¬nite di¬erence

method (12.77), while if φ(t) = t the upwind ¬nite di¬erence scheme (12.78)

(or, equivalently, (12.79)) is obtained. Other choices are admissible as well.

For instance, taking

φ(t) = t ’ 1 + B(2t) (12.83)

where B(t) is the Bernoulli function de¬ned as B(t) = t/(et ’ 1) for t = 0

and B(0) = 1, yields the so called exponential ¬tting ¬nite di¬erence scheme

which is also well known as the Scharfetter-Gummel (SG) method [SG69].

566 12. Two-Point Boundary Value Problems

1.2

1