ńņš. 22 |

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)

4Īµ

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)

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

ńņš. 22 |