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