<< . .

. 80
( : 95)



. . >>

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

<< . .

. 80
( : 95)



. . >>