<<

. 22
( 26)



>>


By Taylor series expansion and recalling (10.66), one obtains

„h (xj ) = ’h’2 [u(xj’1 ) ’ 2u(xj ) + u(xj+1 )] ’ f (xj )
h2 (iv)
= ’u (xj ) ’ f (xj ) + (u (ξj ) + u(iv) (·j )) (12.19)
24
2
h
(u(iv) (ξj ) + u(iv) (·j ))
=
24
for suitable ξj ∈ (xj’1 , xj ) and ·j ∈ (xj , xj+1 ). Upon de¬ning the discrete
maximum norm as
vh h,∞ = max |vh (xj )|,
0¤j¤n

we obtain from (12.19)
f ∞2
¤
„h h (12.20)
h,∞
12
provided that f ∈ C 2 ([0, 1]). In particular, lim „h = 0 and there-
h,∞
h’0
fore the ¬nite di¬erence scheme is consistent with the di¬erential problem
(12.1)-(12.2).

Remark 12.2 Taylor™s expansion of u around xj can also be written as

h2 h3
u(xj ± h) = u(xj ) ± hu (xj ) + u (xj ) ± u (xj ) + R4 (xj ± h)
2 6
538 12. Two-Point Boundary Value Problems

with the following integral form of the remainder
xj +h
(xj + h ’ t)2
R4 (xj + h) = (u (t) ’ u (xj )) dt,
2
xj
xj
(xj ’ h ’ t)2
R4 (xj ’ h) = ’ (u (t) ’ u (xj )) dt.
2
xj ’h

Using the two formulae above, by inspection on (12.18) it is easy to see
that
1
(R4 (xj + h) + R4 (xj ’ h)) .
„h (xj ) = (12.21)
h2
For any integer m ≥ 0, we denote by C m,1 (0, 1) the space of all functions
in C m (0, 1) whose m-th derivative is Lipschitz continuous, i.e.

|v (m) (x) ’ v (m) (y)|
¤ M < ∞.
max
|x ’ y|
x,y∈(0,1),x=y

Looking at (12.21) we see that it su¬ces to assuming that u ∈ C 3,1 (0, 1)
to conclude that
„h h,∞ ¤ M h2
which shows that the ¬nite di¬erence scheme is consistent with the di¬er-
ential problem (12.1)-(12.2) even under a slightly weaker regularity of the
exact solution u.

Remark 12.3 Let e = u ’ uh be the discretization error grid function.
Then,
Lh e = Lh u ’ Lh uh = Lh u ’ fh = „h . (12.22)
It can be shown (see Exercise 5) that

¤3
2 2 2
„h f +f (12.23)
L2 (0,1)
h h

from which it follows that the norm of the discrete second-order derivative
of the discretization error is bounded, provided that the norms of f at the
right-hand side of (12.23) are also bounded.


12.2.2 Convergence Analysis
The ¬nite di¬erence solution uh can be characterized by a discrete Green™s
function as follows. For a given grid point xk de¬ne a grid function Gk ∈ Vh
0

as the solution to the following problem
Lh Gk = ek , (12.24)
12.2 Finite Di¬erence Approximation 539

where ek ∈ Vh satis¬es ek (xj ) = δkj , 1 ¤ j ¤ n ’ 1. It is easy to see
0

that Gk (xj ) = hG(xj , xk ), where G is the Green™s function introduced in
(12.4) (see Exercise 6). For any grid function g ∈ Vh we can de¬ne the grid
0

function
n’1
g(xk )Gk .
wh = Th g, wh = (12.25)
k=1

Then
n’1 n’1
k
Lh wh = g(xk )Lh G = g(xk )ek = g.
k=1 k=1

In particular, the solution uh of (12.10) satis¬es uh = Th f , therefore
n’1 n’1
k
uh = f (xk )G , and uh (xj ) = h G(xj , xk )f (xk ). (12.26)
k=1 k=1

Theorem 12.1 Assume that f ∈ C 2 ([0, 1]). Then, the nodal error e(xj ) =
u(xj ) ’ uh (xj ) satis¬es
h2
u ’ uh h,∞ ¤ f ∞, (12.27)
96
i.e. uh converges to u (in the discrete maximum norm) with second order
with respect to h.
Proof. We start by noticing that, thanks to the representation (12.25), the
following discrete counterpart of (12.5) holds
1
¤
uh f h,∞ . (12.28)
h,∞
8
Indeed, we have
n’1 n’1
|uh (xj )| ¤h G(xj , xk )|f (xk )| ¤ f h G(xj , xk )
h,∞
k=1 k=1
1 1
= f h,∞ xj (1 ’ xj ) ¤ f h,∞
2 8
since, if g = 1, then Th g is such that Th g(xj ) = 1 xj (1 ’ xj ) (see Exercise 7).
2
Inequality (12.28) provides a result of stability in the discrete maximum norm
for the ¬nite di¬erence solution uh . Using (12.22), by the same argument used to
prove (12.28) we obtain
1
¤
e „h h,∞ .
h,∞
8
3
Finally, the thesis (12.27) follows owing to (12.20).

Observe that for the derivation of the convergence result (12.27) we have
used both stability and consistency. In particular, the discretization error
is of the same order (with respect to h) as the consistency error „h .
540 12. Two-Point Boundary Value Problems

12.2.3 Finite Di¬erences for Two-Point Boundary Value
Problems with Variable Coe¬cients
A two-point boundary value problem more general than (12.1)-(12.2) is the
following one
Lu(x) = ’(J(u)(x)) + γ(x)u(x) = f (x) 0 < x < 1,
(12.29)
u(0) = d0 , u(1) = d1
where

J(u)(x) = ±(x)u (x), (12.30)

d0 and d1 are assigned constants and ±, γ and f are given functions that
are continuous in [0, 1]. Finally, γ(x) ≥ 0 in [0, 1] and ±(x) ≥ ±0 > 0 for a
suitable ±0 . The auxiliary variable J(u) is the ¬‚ux associated with u and
very often has a precise physical meaning.
For the approximation, it is convenient to introduce on [0, 1] a new grid
made by the midpoints xj+1/2 = (xj + xj+1 )/2 of the intervals [xj , xj+1 ]
for j = 0, . . . , n ’ 1. Then, a ¬nite di¬erence approximation of (12.29) is
given by: ¬nd uh ∈ Vh such that
Lh uh (xj ) = f (xj ) for all j = 1, . . . , n ’ 1,
(12.31)
uh (x0 ) = d0 , uh (xn ) = d1 ,
where Lh is de¬ned for j = 1, . . . , n ’ 1 as
Jj+1/2 (wh ) ’ Jj’1/2 (wh )
Lh w(xj ) = ’ + γj wj . (12.32)
h
We have de¬ned γj = γ(xj ) and, for j = 0, . . . , n ’ 1, the approximate
¬‚uxes are given by
wj+1 ’ wj
Jj+1/2 (wh ) = ±j+1/2 (12.33)
h
with ±j+1/2 = ±(xj+1/2 ).
The ¬nite di¬erence scheme (12.31)-(12.32) with the approximate ¬‚uxes
(12.33) can still be cast in the form (12.7) by setting

Afd = h’2 tridiagn’1 (a, d, a) + diagn’1 (c) (12.34)

where
T
∈ Rn’2 ,
a = ±1/2 , ±3/2 , . . . , ±n’1/2
T
∈ Rn’1 ,
d = ±1/2 + ±3/2 , . . . , ±n’3/2 + ±n’1/2
T
c = (γ1 , . . . , γn’1 ) ∈ Rn’1 .
12.2 Finite Di¬erence Approximation 541

The matrix (12.34) is symmetric positive de¬nite and is also strictly diag-
onally dominant if γ > 0.
The convergence analysis of the scheme (12.31)-(12.32) can be carried
out by extending straightforwardly the techniques developed in Sections
12.2.1 and 12.2.2.

We conclude this section by addressing boundary conditions that are more
general than those considered in (12.29). For this purpose we assume that

u(0) = d0 , J(u(1)) = g1 ,

where d0 and g1 are two given data. The boundary condition at x = 1 is
called a Neumann condition while the condition at x = 0 (where the value
of u is assigned) is a Dirichlet boundary condition. The ¬nite di¬erence
discretization of the Neumann boundary condition can be performed by
using the mirror imaging technique. For any su¬ciently smooth function
ψ we write its truncated Taylor™s expansion at xn as

h2
ψn’1/2 + ψn+1/2
’ (ψ (·n ) + ψ (ξn ))
ψn =
2 16
for suitable ·n ∈ (xn’1/2 , xn ), ξn ∈ (xn , xn+1/2 ). Taking ψ = J(u) and
neglecting the term containing h2 yields

Jn+1/2 (uh ) = 2g1 ’ Jn’1/2 (uh ). (12.35)

Notice that the point xn+1/2 = xn + h/2 and the corresponding ¬‚ux Jn+1/2
do not really exist (indeed, xn+1/2 is called a “ghost” point), but it is
generated by linear extrapolation of the ¬‚ux at the nearby nodes xn’1/2
and xn . The ¬nite di¬erence equation (12.32) at the node xn reads

Jn’1/2 (uh ) ’ Jn+1/2 (uh )
+ γn un = fn .
h
Using (12.35) to obtain Jn+1/2 (uh ) we ¬nally get the second-order accurate
approximation
±n’1/2
un’1 γn g1 fn
’±n’1/2 + + un = +.
h2 h2 2 h 2
This formula suggests easy modi¬cation of the matrix and right-hand side
entries in the ¬nite di¬erence system (12.7).
For a further generalization of the boundary conditions in (12.29) and
their discretization using ¬nite di¬erences we refer to Exercise 10 where
boundary conditions of the form »u + µu = g at both the endpoints of
(0, 1) are considered for u (Robin boundary conditions).
For a thorough presentation and analysis of ¬nite di¬erence approxima-
tions of two-point boundary value problems, see, e.g., [Str89] and [HGR96].
542 12. Two-Point Boundary Value Problems

12.3 The Spectral Collocation Method
Other discretization schemes can be derived which exhibit the same struc-
ture as the ¬nite di¬erence problem (12.10), with a discrete operator Lh
being de¬ned in a di¬erent manner, though.
Actually, numerical approximations of the second derivative other than
the centred ¬nite di¬erence one can be set up, as described in Section
10.10.3. A noticeable instance is provided by the spectral collocation method.
In that case we assume the di¬erential equation (12.1) to be set on the in-
terval (’1, 1) and choose the nodes {x0 , . . . , xn } to coincide with the n + 1
Legendre-Gauss-Lobatto nodes introduced in Section 10.4. Besides, uh is a
polynomial of degree n. For coherence, we will use throughout the section
the index n instead of h.
The spectral collocation problem reads

¬nd un ∈ P0 : Ln un (xj ) = f (xj ), j = 1, . . . , n ’ 1 (12.36)
n

where P0 is the set of polynomials p ∈ Pn ([0, 1]) such that p(0) = p(1) = 0.
n
Besides, Ln v = LIn v for any continuous function v where In v ∈ Pn is the
interpolant of v at the nodes {x0 , . . . , xn } and L denotes the di¬erential
operator at hand, which, in the case of equation (12.1), coincides with
’d2 /dx2 . Clearly, if v ∈ Pn then Ln v = Lv.
The algebraic form of (12.36) becomes

Asp u = f ,

where uj = un (xj ), fj = f (xj ) j = 1, . . . , n’1 and the spectral collocation
˜ ˜
matrix Asp ∈ R(n’1)—(n’1) is equal to D2 , where D is the matrix obtained
from the pseudo-spectral di¬erentiation matrix (10.73) by eliminating the
¬rst and the n + 1-th rows and columns.
For the analysis of (12.36) we can introduce the following discrete scalar
product
n
(u, v)n = u(xj )v(xj )wj , (12.37)
j=0

where wj are the weights of the Legendre-Gauss-Lobatto quadrature for-
mula (see Section 10.4). Then (12.36) is equivalent to

∀vn ∈ P0 .
(Ln un , vn )n = (f, vn )n (12.38)
n

Since (12.37) is exact for u, v such that uv ∈ P2n’1 (see Section 10.2) then

∀vn ∈ P0 .
2
(Ln vn , vn )n = (Ln vn , vn ) = vn L2 (’1,1) , n

Besides,

(f, vn )n ¤ f ¤
vn 6f vn L2 (’1,1) ,

n n
12.3 The Spectral Collocation Method 543

where f ∞ denotes the maximum of f in [’1, 1] and we have used the

fact that f n ¤ 2 f ∞ and the result of equivalence

vn L2 (’1,1) ¤ vn n ¤ 3 vn L2 (’1,1) , ∀vn ∈ Pn

(see [CHQZ88], p. 286).
Taking vn = un in (12.38) and using the Poincar´ inequality (12.16) we
e
¬nally obtain

un L2 (’1,1) ¤ 6CP f ∞

which ensures that problem (12.36) has a unique solution which is stable.
As for consistency, we can notice that

„n (xj ) = (Ln u ’ f )(xj ) = (’(In u) ’ f )(xj ) = (u ’ In u) (xj )

and this right-hand side tends to zero as n ’ ∞ provided that u ∈
C 2 ([’1, 1]).
Let us now establish a convergence result for the spectral collocation
approximation of (12.1). In the following, C is a constant independent of
n that can assume di¬erent values at di¬erent places.
Moreover, we denote by Hs (a, b), for s ≥ 1, the space of the functions f ∈
C s’1 (a, b) such that f (s’1) is continuous and piecewise di¬erentiable, so
that f (s) exists unless for a ¬nite number of points and belongs to L2 (a, b).
The space Hs (a, b) is known as the Sobolev function space of order s and is
endowed with the norm · Hs (a,b) de¬ned in (10.35).

Theorem 12.2 Let f ∈ Hs (’1, 1) for some s ≥ 1. Then

¤ Cn’s
u ’ un f +u . (12.39)
L2 (’1,1) Hs (’1,1) Hs+1 (’1,1)

Proof. Note that un satis¬es
(un , vn ) = (f, vn )n
1
uvdx is the scalar product of L2 (’1, 1). Similarly, u satis¬es
where (u, v) = ’1

∀v ∈ C 1 ([0, 1]) such that v(0) = v(1) = 0
(u , v ) = (f, v)

(see (12.43) of Section 12.4). Then

((u ’ un ) , vn ) = (f, vn ) ’ (f, vn )n =: E(f, vn ), ∀vn ∈ P0 .
n

It follows that
((u ’ un ) , (u ’ un ) ) = ((u ’ un ) , (u ’ In u) ) + ((u ’ un ) , (In u ’ un ) )
= ((u ’ un ) , (u ’ In u) ) + E(f, In u ’ un ).

We recall the following result (see (10.36))

|E(f, vn )| ¤ Cn’s f vn L2 (’1,1) .
Hs (’1,1)
544 12. Two-Point Boundary Value Problems

Then

|E(f, In u ’ un )| ¤ Cn’s f In u ’ u + u ’ un .
Hs (’1,1) L2 (’1,1) L2 (’1,1)

We recall now the following Young™s inequality (see Exercise 8)
12
ab ¤ µa2 + ∀a, b ∈ R, ∀µ > 0.
b, (12.40)

Using this inequality we obtain
1
(u ’ un ) , (u ’ In u) ¤ (u ’ un ) + (u ’ In u)
2 2
L2 (’1,1) ,
L2 (’1,1)
4
and also (using the Poincar´ inequality (12.16))
e

Cn’s f ¤ C CP n’s f
u ’ un (u ’ un )
Hs (’1,1) Hs (’1,1)
L2 (’1,1) L2 (’1,1)

1
¤ (CCP )2 n’2s f (u ’ un ) 2
2
+ L2 (’1,1) .
Hs (’1,1)
4
Finally,
1 2 ’2s 1
Cn’s f In u ’ u ¤ In u ’ u
2 2
Cn f + L2 (’1,1) .
Hs (’1,1) Hs (’1,1)
L2 (’1,1)
2 2
Using the interpolation error estimate (10.22) for u ’ In u we ¬nally obtain the
3
desired error estimate (12.39).



12.4 The Galerkin Method
We now derive the Galerkin approximation of problem (12.1)-(12.2), which
is the basic ingredient of the ¬nite element method and the spectral method,
widely employed in the numerical approximation of boundary value prob-
lems.


12.4.1 Integral Formulation of Boundary Value Problems
We consider a problem which is slightly more general than (12.1), namely

’(±u ) (x) + (βu )(x) + (γu)(x) = f (x) 0 < x < 1, (12.41)

with u(0) = u(1) = 0, where ±, β and γ are continuous functions on [0, 1]
with ±(x) ≥ ±0 > 0 for any x ∈ [0, 1]. Let us now multiply (12.41) by
a function v ∈ C 1 ([0, 1]), hereafter called a “test function”, and integrate
over the interval [0, 1]
1 1 1 1

f v dx + [±u v]1 ,
±u v dx + βu v dx + γuv dx = 0
0 0 0 0
12.4 The Galerkin Method 545

where we have used integration by parts on the ¬rst integral. If the function
v is required to vanish at x = 0 and x = 1 we obtain
1 1 1 1

±u v dx + βu v dx + γuv dx = f v dx.
0 0 0 0

We will denote by V the test function space. This consists of all functions v
that are continuous, vanish at x = 0 and x = 1 and whose ¬rst derivative is
piecewise continuous, i.e., continuous everywhere except at a ¬nite number
of points in [0, 1] where the left and right limits v’ and v+ exist but do not
necessarily coincide.
V is actually a vector space which is denoted by H1 (0, 1). Precisely,
0

H1 (0, 1) = v ∈ L2 (0, 1) : v ∈ L2 (0, 1), v(0) = v(1) = 0 (12.42)
0

where v is the distributional derivative of v whose de¬nition is given in
Section 12.4.2.

We have therefore shown that if a function u ∈ C 2 ([0, 1]) satis¬es (12.41),
then u is also a solution of the following problem

¬nd u ∈ V : a(u, v) = (f, v) for all v ∈ V, (12.43)
1
f v dx denotes the scalar product of L2 (0, 1) and
where now (f, v) = 0

1 1 1

a(u, v) = ±u v dx + βu v dx + γuv dx (12.44)
0 0 0

is a bilinear form, i.e. it is linear with respect to both arguments u and
v. Problem (12.43) is called the weak formulation of problem (12.1). Since
(12.43) contains only the ¬rst derivative of u it might cover cases in which
a classical solution u ∈ C 2 ([0, 1]) of (12.41) does not exist although the
physical problem is well de¬ned.
If for instance, ± = 1, β = γ = 0, the solution u(x) denotes of the
displacement at point x of an elastic cord having linear density equal to f ,
whose position at rest is u(x) = 0 for all x ∈ [0, 1] and which remains ¬xed
at the endpoints x = 0 and x = 1. Figure 12.2 (right) shows the solution
u(x) corresponding to a function f which is discontinuous (see Figure 12.2,
left). Clearly, u does not exist at the points x = 0.4 and x = 0.6 where f
is discontinuous.
If (12.41) is supplied with non homogeneous boundary conditions, say
u(0) = u0 , u(1) = u1 , we can still obtain a formulation like (12.43) by
proceeding as follows. Let u(x) = xu1 + (1 ’ x)u0 be the straight line that
¯
0 0
interpolates the data at the endpoints, and set u= u(x) ’ u(x). Then u∈ V
¯
546 12. Two-Point Boundary Value Problems

0

’0.005

f(x) ’0.01


u(x)
’0.015

0.6
0.4 1x
0 ’0.02

’0.025

’0.03

’0.035

’1 ’0.04

’0.045

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




FIGURE 12.2. Elastic cord ¬xed at the endpoints and subject to a discontinuous
load f (left). The vertical displacement u is shown on the right

satis¬es the following problem
0 0
¬nd u∈ V : a(u, v) = (f, v) ’ a(¯, v) for all v ∈ V.
u

A similar problem is obtained in the case of Neumann boundary conditions,
say u (0) = u (1) = 0. Proceeding as we did to obtain (12.43), we see
that the solution u of this homogeneous Neumann problem satis¬es the
same problem (12.43) provided the space V is now H1 (0, 1). More general
boundary conditions of mixed type can be considered as well (see Exercise
12).


12.4.2 A Quick Introduction to Distributions
Let X be a Banach space, i.e., a normed and complete vector space. We
say that a functional T : X ’ R is continuous if limx’x0 T (x) = T (x0 )
for all x0 ∈ X and linear if T (x + y) = T (x) + T (y) for any x, y ∈ X and
T (»x) = »T (x) for any x ∈ X and » ∈ R.
Usually, a linear continuous functional is denoted by T, x and the sym-
bol ·, · is called duality. As an example, let X = C 0 ([0, 1]) be endowed
with the maximum norm · ∞ and consider on X the two functionals
de¬ned as
1
T, x = x(0), S, x = x(t) sin(t)dt.
0
It is easy to check that both T and S are linear and continuous functionals
on X. The set of all linear continuous functionals on X identi¬es an abstract
space which is called the dual space of X and is denoted by X .

We then introduce the space C0 (0, 1) (or D(0, 1)) of in¬nitely di¬er-
entiable functions having compact support in [0, 1], i.e., vanishing outside
a bounded open set (a, b) ‚ (0, 1) with 0 < a < b < 1. We say that
vn ∈ D(0, 1) converges to v ∈ D(0, 1) if there exists a closed bounded set
K ‚ (0, 1) such that vn vanishes outside K for each n and for any k ≥ 0
(k)
the derivative vn converges to v (k) uniformly in (0, 1).
12.4 The Galerkin Method 547

The space of linear functionals on D(0, 1) which are continuous with
respect to the convergence introduced above is denoted by D (0, 1) (the
dual space of D(0, 1)) and its elements are called distributions.
We are now in position to introduce the derivative of a distribution. Let
T be a distribution, i.e. an element of D (0, 1). Then, for any k ≥ 0, T (k)
is also a distribution, de¬ned as

∀• ∈ D(0, 1).
T (k) , • = (’1)k T, •(k) , (12.45)

As an example, consider the Heaviside function

x ≥ 0,
1
H(x) =
0 x < 0.

The distributional derivative of H is the Dirac mass δ at the origin, de¬ned
as
v ’ δ(v) = v(0), v ∈ D(R).
From the de¬nition (12.45), it turns out that any distribution is in¬nitely
di¬erentiable; moreover, if T is a di¬erentiable function its distributional
derivative coincides with the usual one.


12.4.3 Formulation and Properties of the Galerkin Method
Unlike the ¬nite di¬erence method which stems directly from the di¬er-
ential (or strong) form (12.41), the Galerkin method is based on the weak
formulation (12.43). If Vh is a ¬nite dimensional vector subspace of V , the
Galerkin method consists of approximating (12.43) by the problem

¬nd uh ∈ Vh : a(uh , vh ) = (f, vh ) ∀vh ∈ Vh . (12.46)

This is a ¬nite dimensional problem. Actually, let {•1 , . . . , •N } denote a
basis of Vh , i.e. a set of N linearly independent functions of Vh . Then we
can write
N
uh (x) = uj •j (x).
j=1


The integer N denotes the dimension of the vector space Vh . Taking vh = •i
in (12.46), it turns out that the Galerkin problem (12.46) is equivalent to
seeking N unknown coe¬cients {u1 , . . . , uN } such that

N
∀i = 1, . . . , N.
uj a(•j , •i ) = (f, •i ) (12.47)
j=1
548 12. Two-Point Boundary Value Problems

We have used the linearity of a(·, ·) with respect to its ¬rst argument, i.e.
N N
a( uj •j , •i ) = uj a(•j , •i ).
j=1 j=1

If we introduce the matrix AG = (aij ), aij = a(•j , •i ) (called the sti¬ness
matrix), the unknown vector u = (u1 , . . . , uN ) and the right-hand side
vector fG = (f1 , . . . , fN ), with fi = (f, •i ), we see that (12.47) is equivalent
to the linear system
AG u = fG . (12.48)
The structure of AG , as well as the degree of accuracy of uh , depends on
the form of the basis functions {•i }, and therefore on the choice of Vh .
We will see two remarkable instances, the ¬nite element method, where
Vh is a space of piecewise polynomials over subintervals of [0, 1] of length not
greater than h which are continuous and vanish at the endpoints x = 0 and
1, and the spectral method in which Vh is a space of algebraic polynomials
still vanishing at the endpoints x = 0, 1.
However, before speci¬cally addressing those cases, we state a couple of
general results that hold for any Galerkin problem (12.46).


12.4.4 Analysis of the Galerkin Method
We endow the space H1 (0, 1) with the following norm
0
±1 1/2
 
|v|H1 (0,1) = |v (x)| dx
2
. (12.49)
 
0

We will address the special case where β = 0 and γ(x) ≥ 0. In the most
general case given by the di¬erential problem (12.41) we shall assume that
the coe¬cients satisfy
1
’ β + γ ≥ 0, ∀x ∈ [0, 1]. (12.50)
2
This ensures that the Galerkin problem (12.46) admits a unique solution
depending continuously on the data. Taking vh = uh in (12.46) we obtain
1 1

±0 |uh |2 1 (0,1) ¤ γuh uh dx = (f, uh ) ¤ f
±uh uh dx + uh L2 (0,1) ,
L2 (0,1)
H
0 0

where we have used the Cauchy-Schwarz inequality (8.29) to set the right-
hand side inequality. Owing to the Poincar´ inequality (12.16) we conclude
e
that
CP
|uh |H1 (0,1) ¤ f L2 (0,1) . (12.51)
±0
12.4 The Galerkin Method 549

Thus, the norm of the Galerkin solution remains bounded (uniformly with
respect to the dimension of the subspace Vh ) provided that f ∈ L2 (0, 1).
Inequality (12.51) therefore represents a stability result for the solution of
the Galerkin problem.
As for convergence, we can prove the following result.

’1 2
Theorem 12.3 Let C = ±0 ( ± + CP γ ∞ ); then, we have



|u ’ uh |H1 (0,1) ¤ C min |u ’ wh |H1 (0,1) . (12.52)
wh ∈Vh

Proof. Subtracting (12.46) from (12.43) (where we use vh ∈ Vh ‚ V ), owing to
the bilinearity of the form a(·, ·) we obtain

a(u ’ uh , vh ) = 0 ∀vh ∈ Vh . (12.53)

Then, setting e(x) = u(x) ’ uh (x), we deduce

±0 |e|2 1 (0,1) ¤ a(e, e) = a(e, u ’ wh ) + a(e, wh ’ uh ) ∀wh ∈ Vh .
H


The last term is null due to (12.53). On the other hand, still by the Cauchy-
Schwarz inequality we obtain
1 1

a(e, u ’ wh ) ±e (u ’ wh ) dx + γe(u ’ wh ) dx
=
0 0
¤± ∞ |e|H1 (0,1) |u ’ wh |H1 (0,1) + γ u ’ wh
e L2 (0,1) .
∞ L2 (0,1)


The desired result (12.52) now follows by using again the Poincar´ inequality for
e
both e L2 (0,1) and u ’ wh L2 (0,1) . 3

The previous results can be obtained under more general hypotheses on
problems (12.43) and (12.46). Precisely, we can assume that V is a Hilbert
space, endowed with norm · V , and that the bilinear form a : V — V ’ R
satis¬es the following properties:

∃±0 > 0 : a(v, v) ≥ ±0 v ∀v ∈ V (coercivity),
2
(12.54)
V



∃M > 0 : |a(u, v)| ¤ M u ∀u, v ∈ V (continuity).
v (12.55)
V V


Moreover, the right hand side (f, v) satis¬es the following inequality

|(f, v)| ¤ K v ∀v ∈ V.
V


Then both problems (12.43) and (12.46) admit unique solutions that satisfy

K K
¤ ¤
u , uh .
V V
±0 ±0
550 12. Two-Point Boundary Value Problems

This is a celebrated result which is known as the Lax-Milgram Lemma (for
its proof see, e.g., [QV94]). Besides, the following error inequality holds
M
u ’ uh ¤ min u ’ wh . (12.56)
V V
±0 wh ∈Vh
The proof of this last result, which is known as C´a™s Lemma, is very similar
e
to that of (12.52) and is left to the reader.

We now wish to notice that, under the assumption (12.54), the matrix
introduced in (12.48) is positive de¬nite. To show this, we must check that
vT Bv ≥ 0 ∀v ∈ RN and that vT Bv = 0 ” v = 0 (see Section 1.12).
Let us associate with a generic vector v = (vi ) of RN the function vh =
N
j=1 vj •j ∈ Vh . Since the form a(·, ·) is bilinear and coercive we get

N N N N
T
v AG v = vi aij vj = vi a(•j , •i )vj
« 
j=1 i=1 j=1 i=1
N N N N
a(vj •j , vi •i ) = a  vi •i 
= vj •j ,
j=1 i=1 j=1 i=1
= a(vh , vh ) ≥ ± ≥ 0.
vh 2
V

Moreover, if vT AG v = 0 then also vh 2
= 0 which implies vh = 0 and
V
thus v = 0.

It is also easy to check that the matrix AG is symmetric i¬ the bilinear
form a(·, ·) is symmetric.
For example, in the case of problem (12.41) with β = γ = 0 the ma-
trix AG is symmetric and positive de¬nite (s.p.d.) while if β and γ are
nonvanishing, AG is positive de¬nite only under the assumption (12.50).
If AG is s.p.d. the numerical solution of the linear system (12.48) can be
e¬ciently carried out using direct methods like the Cholesky factorization
(see Section 3.4.2) as well as iterative methods like the conjugate gradient
method (see Section 4.3.4). This is of particular interest in the solution of
boundary value problems in more than one space dimension (see Section
12.6).


12.4.5 The Finite Element Method
The ¬nite element method (FEM) is a special technique for constructing
a subspace Vh in (12.46) based on the piecewise polynomial interpolation
considered in Section 8.3. With this aim, we introduce a partition Th of
[0,1] into n subintervals Ij = [xj , xj+1 ], n ≥ 2, of width hj = xj+1 ’ xj ,
j = 0, . . . , n ’ 1, with
0 = x0 < x1 < . . . < xn’1 < xn = 1
12.4 The Galerkin Method 551

and let h = max(hj ). Since functions in H1 (0, 1) are continuous it makes
0
Th
sense to consider for k ≥ 1 the family of piecewise polynomials Xh intro-
k

duced in (8.22) (where now [a, b] must be replaced by [0, 1]). Any function
vh ∈ Xh is a continuous piecewise polynomial over [0, 1] and its restriction
k

over each interval Ij ∈ Th is a polynomial of degree ¤ k. In the following
we shall mainly deal with the cases k = 1 and k = 2.
Then, we set
k,0
Vh = Xh = vh ∈ Xh : vh (0) = vh (1) = 0 .
k
(12.57)

The dimension N of the ¬nite element space Vh is equal to nk ’ 1. In the
following the two cases k = 1 and k = 2 will be examined.
To assess the accuracy of the Galerkin FEM we ¬rst notice that, thanks
to C´a™s lemma (12.56), we have
e

min u ’ wh ¤ u ’ Πk u (12.58)
H1 (0,1) H1 (0,1)
h
0 0
wh ∈Vh


where Πk u is the interpolant of the exact solution u ∈ V of (12.43) (see
h
Section 8.3). From inequality (12.58) we conclude that the matter of esti-
mating the Galerkin approximation error u ’ uh H1 (0,1) is turned into the
0
estimate of the interpolation error u ’ Πk u H1 (0,1) . When k = 1, using
h 0
(12.56) and (8.27) we obtain

M
u ’ uh ¤ Ch u
H1 (0,1) H2 (0,1)
0
±0

provided that u ∈ H2 (0, 1). This estimate can be extended to the case k > 1
as stated in the following convergence result (for its proof we refer, e.g., to
[QV94], Theorem 6.2.1).

Property 12.1 Let u ∈ H1 (0, 1) be the exact solution of (12.43) and uh ∈
0
Vh its ¬nite element approximation using continuous piecewise polynomials
of degree k ≥ 1. Assume also that u ∈ Hs (0, 1) for some s ≥ 2. Then the
following error estimate holds
M
u ’ uh ¤ Chl u (12.59)
H1 (0,1) Hl+1 (0,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 (0,1) . (12.60)
L2 (0,1)



The estimate (12.59) shows that the Galerkin method is convergent, i.e. the
approximation error tends to zero as h ’ 0 and the order of convergence is
552 12. Two-Point Boundary Value Problems

k. We also see that there is no convenience in increasing the degree k of the
¬nite element approximation if the solution u is not su¬ciently smooth. In
this respect l is called a regularity threshold. The obvious alternative to gain
accuracy in such a case is to reduce the stepzise h. Spectral methods, which
will be considered in Section 12.4.7, instead pursue the opposite strategy
(i.e. increasing the degree k) and are thus ideally suited to approximating
problems with highly smooth solutions.
An interesting situation is that where the exact solution u has the min-
imum regularity (s = 1). In such a case, C´a™s lemma ensures that the
e
Galerkin FEM is still convergent since as h ’ 0 the subspace Vh becomes
dense into V . However, the estimate (12.59) is no longer valid so that
it is not possible to establish the order of convergence of the numerical
method. Table 12.1 summarizes the orders of convergence of the FEM for
k = 1, . . . , 4 and s = 1, . . . , 5.

k s=1 s=2 s=3 s=4 s=5
h1 h1 h1 h1
1 only convergence
h2
h1 h2 h2
2 only convergence
h3
h1 h2 h3
3 only convergence
h4
h1 h2 h3
4 only convergence
TABLE 12.1. Order of convergence of the FEM as a function of k (the degree of
interpolation) and s (the Sobolev regularity of the solution u)



Let us now focus on how to generate a suitable basis {•j } for the ¬nite
k
element space Xh in the special cases k = 1 and k = 2. The basic point is
to choose appropriately a set of degrees of freedom for each element Ij of
the partition Th (i.e., the parameters which permit uniquely identifying a
k k
function in Xh ). The generic function vh in Xh can therefore be written as
nk
vh (x) = vi •i (x)
i=0

where {vi } denote the set of the degrees of freedom of vh and the basis
functions •i (which are also called shape functions) are assumed to satisfy
the Lagrange interpolation property •i (xj ) = δij , i, j = 0, . . . , n, where δij
is the Kronecker symbol.

1
The space Xh
This space consists of all continuous and piecewise linear functions over
the partition Th . Since a unique straight line passes through two distinct
nodes the number of degrees of freedom for vh is equal to the number
n + 1 of nodes in the partition. As a consequence, n + 1 shape functions
12.4 The Galerkin Method 553

1
•i , i = 0, . . . , n, are needed to completely span the space Xh . The most
natural choice for •i , i = 1, . . . , n ’ 1, is
±
 x ’ xi’1 for xi’1 ¤ x ¤ xi ,

 x ’x
i

 i’1

xi+1 ’ x
•i (x) = (12.61)
for xi ¤ x ¤ xi+1 ,
x
 i+1 ’ xi





0 elsewhere.

The shape function •i is thus piecewise linear over Th , its value is 1 at
the node xi and 0 at all the other nodes of the partition. Its support (i.e.,
the subset of [0, 1] where •i is nonvanishing) consists of the union of the
intervals Ii’1 and Ii if 1 ¤ i ¤ n ’ 1 while it coincides with the interval I0
(respectively In’1 ) if i = 0 (resp., i = n). The plots of •i , •0 and •n are
shown in Figure 12.3.




1
•0 •i •n




x0 x 1 xi’1 xi xi+1 xn’1 xn = 1
1
FIGURE 12.3. Shape functions of Xh associated with internal and boundary
nodes


For any interval Ii = [xi , xi+1 ], i = 0, . . . , n ’ 1, the two basis functions •i
and •i+1 can be regarded as the images of two “reference” shape functions
•0 and •1 (de¬ned over the reference interval [0, 1]) through the linear
a¬ne mapping φ : [0, 1] ’ Ii

x = φ(ξ) = xi + ξ(xi+1 ’ xi ), i = 0, . . . , n ’ 1. (12.62)

De¬ning •0 (ξ) = 1 ’ ξ, •1 (ξ) = ξ, the two shape functions •i and •i+1
can be constructed over the interval Ii as

•i (x) = •0 (ξ(x)), •i+1 (x) = •1 (ξ(x))

where ξ(x) = (x ’ xi )/(xi+1 ’ xi ) (see Figure 12.4).
554 12. Two-Point Boundary Value Problems




1 1
φ
’’
•1 •i+1



x
ξ xi xi+1
0 1




FIGURE 12.4. Linear a¬ne mapping φ from the reference interval to the generic
interval of the partition

2
The space Xh
The generic function vh ∈ Xh is a piecewise polynomial of degree 2 over
2

each interval Ii . As such, it can be uniquely determined once three values
of it at three distinct points of Ii are assigned. To ensure continuity of
vh over [0, 1] the degrees of freedom are chosen as the function values at
the nodes xi of Th , i = 0, . . . , n, and at the midpoints of each interval Ii ,
i = 0, . . . , n ’ 1, for a total number equal to 2n + 1. It is convenient to
label the degrees of freedom and the corresponding nodes in the partition
starting from x0 = 0 until x2n = 1 in such a way that the midpoints of
each interval correspond to the nodes with odd index while the endpoints
of each interval correspond to the nodes with even index.
The explicit expression of the single shape function is
±
 (x ’ xi’1 )(x ’ xi’2 ) for xi’2 ¤ x ¤ xi ,

 (x ’ x )(x ’ x )


 i i’1 i i’2

(xi+1 ’ x)(xi+2 ’ x)
(i even) •i (x) = (12.63)
 (xi+1 ’ xi )(xi+2 ’ xi ) for xi ¤ x ¤ xi+2 ,






0 elsewhere,

±
 (xi+1 ’ x)(x ’ xi’1 )
 for xi’1 ¤ x ¤ xi+1 ,
(xi+1 ’ xi )(xi ’ xi’1 )
(i odd) •i (x) = (12.64)


0 elsewhere.
Each basis function enjoys the property that •( xj ) = δij , i, j = 0, . . . , 2n.
2
The shape functions for Xh on the reference interval [0, 1] are
•0 (ξ) = (1 ’ ξ)(1 ’ 2ξ), •1 (ξ) = 4(1 ’ ξ)ξ, •2 (ξ) = ξ(2ξ ’ 1) (12.65)
12.4 The Galerkin Method 555

and are shown in Figure 12.5. As in the case of piecewise linear ¬nite
1
elements of Xh the shape functions (12.63) and (12.64) are the images of
(12.65) through the a¬ne mapping (12.62). Notice that the support of the
basis function •2i+1 associated with the midpoint x2i+1 coincides with the
interval to which the midpoint belongs. Due to its shape •2i+1 is usually
referred to as bubble function.




•1
•0
•2




ξ
1
0 0.5
2
FIGURE 12.5. Basis functions of Xh on the reference interval

So far, we have considered only lagrangian-type shape functions. If this
constraint is removed other kind of bases can be derived. A notable example
(on the reference interval) is given by

ψ0 (ξ) = 1 ’ ξ, ψ1 (ξ) = (1 ’ ξ)ξ, ψ2 (ξ) = ξ. (12.66)

This basis is called hierarchical since it is generated using the shape func-
2
tions of the subspace having an immediately lower dimension than Xh (i.e.
Xh ). Precisely, the bubble function ψ1 ∈ Xh is added to the shape func-
1 2
1
tions ψ0 and ψ2 which belong to Xh . Hierarchical bases can be of some
interest in numerical computations if the local degree of the interpolation
is adaptively increased (p-type adaptivity).
2
To check that (12.66) forms a basis for Xh we must verify that its func-
tions are linearly independent, i.e.

∀ξ ∈ [0, 1] ”
±0 ψ0 (ξ) + ±1 ψ1 (ξ) + ±2 ψ2 (ξ) = 0, ±0 = ±1 = ±2 = 0.

In our case this holds true since if
2
±i ψi (ξ) = ±0 + ξ(±1 ’ ±0 + ±2 ) ’ ±1 ξ 2 = 0, ∀ξ ∈ [0, 1]
i=0

then necessarily ±0 = 0, ±1 = 0 and thus ±2 = 0.
556 12. Two-Point Boundary Value Problems

A procedure analogous to that examined in the sections above can be used
k
in principle to construct a basis for every subspace Xh with k being arbi-
trary. However, it is important to remember that an increase in the degree k
of the polynomial approximation gives rise to an increase of the number of
degrees of freedom of the FEM and, as a consequence, of the computational
cost required for the solution of the linear system (12.48).

Let us now examine the structure and the basic properties of the sti¬ness
matrix associated with system (12.48) in the case of the ¬nite element
method (AG = Afe ).
k
Since the ¬nite element basis functions for Xh have a local support, Afe is
sparse. In the particular case k = 1, the support of the shape function •i is
the union of the intervals Ii’1 and Ii if 1 ¤ i ¤ n ’ 1, and it coincides with
the interval I0 (respectively In’1 ) if i = 0 (resp., i = n). As a consequence,
for a ¬xed i = 1, . . . , n ’ 1, only the shape functions •i’1 and •i+1 have a
nonvanishing support intersection with that of •i , which implies that Afe
is tridiagonal since aij = 0 if j ∈ {i ’ 1, i, i + 1}. In the case k = 2 one
concludes with an analogous argument that Afe is a pentadiagonal matrix.
The condition number of Afe is a function of the grid size h; indeed,

A’1 = O(h’2 )
K2 (Afe ) = Afe 2 2
fe

(for the proof, see [QV94], Section 6.3.2), which demonstrates that the
conditioning of the ¬nite element system (12.48) grows rapidly as h ’
0. This is clearly con¬‚icting with the need of increasing the accuracy of
the approximation and, in multidimensional problems, demands suitable
preconditioning techniques if iterative solvers are used (see Section 4.3.2).

Remark 12.4 (Elliptic problems of higher order) The Galerkin me-
thod in general, and the ¬nite element method in particular, can also be
applied to other type of elliptic equations, for instance to those of fourth
order. In that case, the numerical solution (as well as the test functions)
should be continuous together with their ¬rst derivative. An example has
been illustrated in Section 8.8.1.


12.4.6 Implementation Issues
In this section we implement the ¬nite element (FE) approximation with
piecewise linear elements (k = 1) of the boundary value problem (12.41)
(shortly, BVP) with non homogeneous Dirichlet boundary conditions.
Here is the list of the input parameters of Program 94: Nx is the number
of grid subintervals; I is the interval [a, b], alpha, beta, gamma and f are
the macros corresponding to the coe¬cients in the equation, bc=[ua,ub]
is a vector containing the Dirichlet boundary conditions for u at x = a and
x = b and stabfun is an optional string variable. It can assume di¬erent
12.4 The Galerkin Method 557

values, allowing the user to select the desired type of arti¬cial viscosity that
may be needed for dealing with the problems addressed in Section 12.5.

Program 94 - ellfem : Linear FE for two-point BVPs

function [uh,x] = ellfem(Nx,I,alpha,beta,gamma,f,bc,stabfun)
a = I(1); b = I(2); h = (b-a)/Nx; x = [a+h/2:h:b-h/2];
alpha = eval(alpha); beta = eval(beta); gamma = eval(gamma);
f = eval(f); rhs = 0.5*h*(f(1:Nx-1)+f(2:Nx));
if nargin == 8
[Afe,rhsbc] = femmatr(Nx,h,alpha,beta,gamma,stabfun);
else
[Afe,rhsbc] = femmatr(Nx,h,alpha,beta,gamma);
end
[L,U,P] = lu(Afe);
rhs(1) = rhs(1)-bc(1)*(-alpha(1)/h-beta(1)/2+h*gamma(1)/3+rhsbc(1));
rhs(Nx-1) = rhs(Nx-1)-bc(2)*(-alpha(Nx)/h+beta(Nx)/2+h*gamma(Nx)/3+rhsbc(2));
rhs = P*rhs™; z = L \ rhs;
w = U \ z;
uh = [bc(1), w™, bc(2)]; x = [a:h:b];

Program 95 computes the sti¬ness matrix Afe ; with this aim, the coe¬cients
±, β and γ and the forcing term f are replaced by piecewise constant
functions on each mesh subinterval and the remaining integrals in (12.41),
involving the basis functions and their derivatives, are evaluated exactly.

Program 95 - femmatr : Construction of the sti¬ness matrix

function [Afe,rhsbc] = femmatr(Nx,h,alpha,beta,gamma,stabfun)
for i=2:Nx
dd(i-1)=(alpha(i-1)+alpha(i))/h; dc(i-1)=-(beta(i)-beta(i-1))/2;
dr(i-1)=h*(gamma(i-1)+gamma(i))/3;
if i > 2
ld(i-2) = -alpha(i-1)/h; lc(i-2)=-beta(i-1)/2;
lr(i-2) = h*gamma(i-1)/6;
end
if i < Nx
ud(i-1) = -alpha(i)/h; uc(i-1)=beta(i)/2;
ur(i-1) = h*gamma(i)/6;
end
end
Kd=spdiags([[ld 0]™,dd™,[0 ud]™],-1:1,Nx-1,Nx-1);
Kc=spdiags([[lc 0]™,dc™,[0 uc]™],-1:1,Nx-1,Nx-1);
Kr=spdiags([[lr 0]™,dr™,[0 ur]™],-1:1,Nx-1,Nx-1);
Afe=Kd+Kc+Kr;
if nargin == 6
s=[™[Ks,rhsbc]=™,stabfun,™(Nx,h,alpha,beta);™]; eval(s)
Afe = Afe + Ks;
else
558 12. Two-Point Boundary Value Problems

rhsbc = [0, 0];
end

The H1 -norm of the error can be computed by calling Program 96, which
must be supplied by the macros u and ux containing the expression of the
exact solution u and of u . The computed numerical solution is stored in
the output vector uh, while the vector coord contains the grid coordinates
and h is the mesh size. The integrals involved in the computation of the
H1 -norm of the error are evaluated using the composite Simpson formula
(9.17).

Program 96 - H1error : Computation of the H1 -norm of the error
function [L2err,H1err]=H1error(coord,h,uh,u,udx)
nvert=max(size(coord)); x=[]; k=0;
for i = 1:nvert-1
xm = (coord(i+1)+coord(i))*0.5;
x = [x, coord(i),xm]; k = k + 2;
end
ndof = k+1; x (ndof) = coord (nvert);
uq = eval(u); uxq = eval(udx); L2err = 0; H1err = 0;
for i=1:nvert-1
L2err = L2err + (h/6)*((uh(i)-uq(2*i-1))ˆ2+...
4*(0.5*uh(i)+0.5*uh(i+1)-uq(2*i))ˆ2+(uh(i+1)-uq(2*i+1))ˆ2);
H1err = H1err + (1/(6*h))*((uh(i+1)-uh(i)-h*uxq(2*i-1))ˆ2+...
4*(uh(i+1)-uh(i)-h*uxq(2*i))ˆ2+(uh(i+1)-uh(i)-h*uxq(2*i+1))ˆ2);
end
H1err = sqrt(H1err + L2err); L2err = sqrt(L2err);


Example 12.1 We assess the accuracy of the ¬nite element solution of the fol-
lowing problem. Consider a thin rod of length L whose temperature at x = 0 is
¬xed to t0 while the other endpoint x = L is thermally isolated. Assume that the
rod has a cross-section with constant area equal to A and that the perimeter of
A is p.
The temperature u of the rod at a generic point x ∈ (0, L) is governed by the
following boundary value problem with mixed Dirichlet-Neumann conditions
±
 ’µAu + σpu = 0 x ∈ (0, L),
(12.67)
 u(0) = u , u (L) = 0,
0


where µ denotes the thermal conductivity and σ is the convective transfer coef-
¬cient. The exact solution of the problem is the (smooth) function

cosh[m(L ’ x)]
u(x) = u0 ,
cosh(mL)

where m = σp/µA. We solve the problem by using linear and quadratic ¬nite
elements (k = 1 and k = 2) on a grid with uniform size. In the numerical
12.4 The Galerkin Method 559

computations we assume that the length of the rod is L = 100cm and that the
rod has a circular cross-section of radius 2cm (and thus, A = 4πcm2 , p = 4πcm).
We also set u0 = 10—¦ C, σ = 2 and µ = 200.
Figure 12.6 (left) shows the behavior of the error in the L2 and H1 norms for
the linear and quadratic elements, respectively. Notice the excellent agreement
between the numerical results and the expected theoretical estimates (12.59) and
(12.60), i.e., the orders of convergence in the L2 norm and the H1 norm tend
respectively to k + 1 and k if ¬nite elements of degree k are employed, since the

exact solution is smooth.

1
10
2
10
0
10
0
10
’1
10 ’2
10

’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

<<

. 22
( 26)



>>