<< . .

. 81
( : 95)



. . >>

’2 ’4
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)

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 =

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


<< . .

. 81
( : 95)



. . >>