¬ ·

¬ ·

A = ¬ am+1,1 am+1,2 · · · am+1,m+1 . . . .. .. · ∈ RM—M

. .

¬ ·

¬ ·

¬ ·

.. .. .. .. .. ..

¬ ·

. . . . . .

¬ ·

¬ ·

¬ ·

.. .. .. .. ..

. . . . .

0

aM,M’m · · · aM,M’1 aM,M

is mapped to the matrix

«

··· ···

0 0 a11

¬ ·

···

0 0 a21 a22

¬ ·

¬ ·

. . .

¬ ·

. . .

. . .

¬ ·

¬ ·

··· ···

0 am,1 am,m

¬ ·

· ∈ RM—(m+1) .

B = ¬ am+1,1 ··· · · · am+1,m am+1,m+1

¬ ·

¬ ·

. . . . .

¬ ·

. . . . .

. . . . .

¬ ·

¬ ·

. . . . .

. . . . .

. . . . .

aM,M’m · · · · · · aM,M’1 aM,M

The unused elements of B, i.e., (B)ij for i = 1, . . . , m, j = 1, . . . , m + 1 ’ i,

are here ¬lled with zeros.

For a general band matrix, the matrix B ∈ RM—(2m+1) obtained by the

above conversion has the following form:

«

··· ···

0 0 a11 a12 a1,m+1

¬ ·

··· ··· ···

0 a21 a22 a2,m+2

¬ ·

¬ ·

. . . . . .

¬ ·

. . . . . .

¬ ·

. . . . . .

¬ ·

··· ··· ··· ···

¬ ·

0 am,1 am,2m

¬ ·

··· ··· ··· ··· ···

¬ ·

am+1,1 am+1,2m+1

B =¬ ·.

¬ ·

. . . . . . .

. . . . . . .

¬ ·

. . . . . . .

¬ ·

¬ aM ’m,M ’2m ·

··· ··· ··· ··· ··· aM ’m,M

¬ ·

¬ aM ’m+1,M ’2m+1 · · · ·

··· ··· · · · aM ’m+1,M 0

¬ ·

¬ ·

. . . . . .

. . . . . .

. . . . . .

··· ··· ···

aM,M ’m aM,M 0 0

Here, in the right lower part of the matrix, a further sector of unused

elements arose, which is also ¬lled with zeros.

If the storage is based on the hull, additionally a pointer ¬eld is needed,

which points to the diagonal elements, for example. If the matrix is sym-

2.5. Direct Methods for Sparse Linear Systems 87

metric, again the storage of the lower triangular matrix is su¬cient. For

the matrix A from Figure 2.16 under the assumption that A is symmetric,

the pointer ¬eld could act as shown in Figure 2.17.

a11 a22 a31 a32 a33 a43 a44 a52 a53 a54

a55

$

X

$$$

rrrr T

$$

rr rr $$

$$

rr rr $$

1 2 5 7 11

Figure 2.17. Linear storage of the hull.

Coupled Assembling and Decomposition

A formerly popular method, the so-called frontal method, performs

simultaneously assembling and the Cholesky factorization.

We consider this method for the example of the sti¬ness matrix Ah =

(aij ) ∈ RM—M with bandwidth m (with the original numbering).

The method is based on the kth step of the Gaussian or Cholesky method

(cf. Figure 2.18).

k

0

k

Bk

0

Figure 2.18. kth step of the Cholesky method.

Only the entries of Bk are to be changed, i.e., only those elements aij

with k ¤ i, j ¤ k + m. The corresponding formula is

(k)

aik

(k+1) (k) (k)

’

aij = aij akj , i, j = k + 1, . . . , k + m . (2.58)

(k)

akk

Here, the upper indices indicate the steps of the elimination method, which

we store in aij . The entries aij are generated by summation of entries of

88 2. Finite Element Method for Poisson Equation

the element sti¬ness matrix of those elements K that contain nodes with

the indices i, j.

(k) (k)

Furthermore, to perform the elimination step (2.58), only aik , akj for

(k)

i, j = k, . . . , k + m must be completely assembled; aij , i, j = k + 1, . . . , k +

(k) (k+1) (k+1) (k+1)

m, can be replaced by aij if aij

˜ is later de¬ned by aij := aij

˜ +

(k) (k)

aij ’ aij . That is, for the present, aij needs to consist of only a few

˜

contributions of elements K with nodes i, j in K.

From these observations, the following algorithm is obtained. The kth

step for k = 1, . . . , M reads as follows:

• Assemble all of the missing contributions of elements K that contain

the node with index k.

• Compute A(k+1) by modi¬cation of the entries of Bk according to

(2.58).

• Store the kth row of A(k+1) , also out of the main memory.

• De¬ne Bk+1 (by a south-east shift).

Here the assembling is node-based and not element-based.

The advantage of this method is that Ah need not be completely assem-

bled and stored in the main memory, but only a matrix Bk ∈ R(m+1)—(m+1) .

Of course, if M is not too large, there may be no advantage.

Bandwidth Reduction

The complexity, i.e., the number of operations, is crucial for the application

of a particular method:

The Cholesky method, applied to a symmetric matrix A ∈ RM—M with

bandwidth m, requires O(m2 M ) operations in order to compute L.

However, the bandwidth m of the sti¬ness matrix depends on the num-

bering of the nodes. Therefore, a numbering is to be found where the

number m is as small as possible.

We want to consider this again for the example of the Poisson equation on

the rectangle with the discretization according to Figure 2.9. Let the inte-

rior nodes have the coordinates (ih, jh) with i = 1, . . . k ’1, j = 1, . . . , l’1.

The discretization corresponds to the ¬nite di¬erence method introduced

beginning with (1.10); i.e., the bandwidth is equal to k ’ 1 for a rowwise

numbering or l ’ 1 for a columnwise numbering.

For k l or k l, this fact results in a large di¬erence of the bandwidth

m or of the pro¬le (of the left triangle), which is of size (k ’ 1)(l ’ 1)(m + 1)

except for a term of m2 . Therefore, the columnwise numbering is preferred

for k l; the rowwise numbering is preferred for k l.

For a general domain „¦, a numbering algorithm based on a given tri-

angulation Th and on a basis {•i } of Vh is necessary with the following

properties:

2.5. Direct Methods for Sparse Linear Systems 89

The structure of A resulting from the numbering must be such that the

band or the pro¬le of A is as small as possible. Furthermore, the numbering

algorithm should yield the numbers m(A) or fi (A), mi (A) such that the

matrix A can also be assembled using the element matrices A(k) .

Given a triangulation Th and a corresponding basis •i 1 ¤ i ¤ M of

Vh , we start with the assignment of some graph G to this triangulation as

follows:

The nodes of G coincide with the nodes {a1 , . . . , aM } of the triangulation.

The de¬nition of its edges is:

⇐’

(ai , aj ) is an edge of G

there exists a K ∈ Th such that •i |K ≡ 0 , •j |K ≡ 0 .

In Figure 2.19 some examples are given, where the example (2) will be

introduced in Section 3.3.

triangulation:

. . . . . .

(1) (2)

. . . . .

linear ansatz on triangle (bi)linear ansatz on quadrilateral

Graph:

. . . . . .

. . . . .

Figure 2.19. Triangulation and assigned graph.

If several degrees of freedom are assigned to some node of the triangu-

lation Th , then also in G several nodes are assigned to it. This is the case,

for example, if so-called Hermite elements are considered, which will be

introduced in Section 3.3. The costs of administration are small if the same

number of degrees of freedom is assigned to all nodes of the triangulation.

An often-used numbering algorithm is the Cuthill“McKee method. This

algorithm operates on the graph G just de¬ned. Two nodes ai , aj of G are

called neighboured if (ai , aj ) is an edge of G. The degree of a node ai of G

is de¬ned as the number of neighbours of ai .

The kth step of the algorithm for k = 1, . . . , M has the following form:

k = 1: Choose a starting node, which gets the number 1. This starting

node forms the level 1.

k > 1: If all nodes are already numbered, the algorithm is terminated.

Otherwise, the level k is formed by taking all the nodes that are not num-

90 2. Finite Element Method for Poisson Equation

bered yet and that are neighbours of a node of level k ’ 1. The nodes of

level k will be consecutively numbered.

Within a level, we can sort, for example, by the degree, where the node

with the smallest degree is numbered ¬rst.

The reverse Cuthill“McKee method consists of the above method and the

inversion of the numbering at the end; i.e.,

new node number = M + 1 ’ old node number .

This corresponds to a re¬‚ection of the matrix at the counterdiagonal. The

bandwidth does not change by the inversion, but the pro¬le may diminish

drastically for many examples (cf. Figure 2.20).

* *** *** *

* *** *

* *** *

* *** *

* *** *

* *** *** *

Figure 2.20. Change of the hull by re¬‚ection at the counterdiagonal.

The following estimate holds for the bandwidth m of the numbering

created by the Cuthill“McKee algorithm:

D+i

¤ m ¤ max (Nk’1 + Nk ’ 1) .

2 2¤k¤ν

Here D is the maximum degree of a node of G, ν is the number of levels, and

Nk is the number of nodes of level k. The number i is equal to 0 if D is even,

and i is equal to 1 if D is odd. The left-hand side of the above inequality is

easy to understand by means of the following argument: To reach a minimal

bandwidth, all nodes that are neighbours of ai in the graph G should also

be neighbours of ai in the numbering. Then the best situation is given if the

neighboured nodes would appear uniformly immediately before and after

ai . If D is odd, then one side has one node more than the other.

To verify the right-hand side, consider a node ai that belongs to level

k ’ 1 as well as a node aj that is a neighbour of ai in the graph G and that

is not yet numbered in level k ’ 1. Therefore, aj will get a number in the

kth step. The largest bandwidth is obtained if ai is the ¬rst node of the

numbering of level k ’ 1 and if aj is the last node of level k. Hence exactly

(Nk’1 ’ 1) + (Nk ’ 1) nodes lie between both of these; i.e., their distance

in the numbering is Nk’1 + Nk ’ 1.

It is favourable if the number ν of levels is as large as possible and if all

the numbers Nk are of the same size, if possible. Therefore, the starting

node should be chosen “at one end” of the graph G if possible; if all the

2.5. Direct Methods for Sparse Linear Systems 91

˜ ˜

starting nodes are to be checked, the expense will be O(M M ), where M

is the number of edges of G. One possibility consists in choosing a node

with minimum degree for the starting node. Another possibility is to let

the algorithm run once and then to choose the last-numbered node as the

starting node.

If a numbering is created by the (reverse) Cuthill“McKee algorithm, we

can try to improve it “locally”, i.e., by exchanging particular nodes.

Exercise

2.8 Show that the number of arithmetic operations for the Cholesky

method for an M — M matrix with bandwidth m has order M m2 /2;

additionally, M square roots have to be calculated.

3

The Finite Element Method for Linear

Elliptic Boundary Value Problems of

Second Order

3.1 Variational Equations and Sobolev Spaces

We now continue the de¬nition and analysis of the “correct” function spaces

that we began in (2.17)“(2.20). An essential assumption ensuring the exis-

tence of a solution of the variational equation (2.13) is the completeness of

the basic space (V, · ). In the concrete case of the Poisson equation the

“preliminary” function space V according to (2.7) can be equipped with

the norm · 1 , de¬ned in (2.19), which has been shown to be equivalent to

the norm · a , given in (2.6) (see (2.46)). If we consider the minimization

problem (2.14), which is equivalent to the variational equation, the func-

tional F is bounded from below such that the in¬mum assumes a ¬nite

value and there exists a minimal sequence (vn )n in V , that is, a sequence

with the property

lim F (vn ) = inf F (v) v ∈ V .

n’∞

The form of F also implies that (vn )n is a Cauchy sequence. If this sequence

converges to an element v ∈ V , then, due to the continuity of F with respect

to · , it follows that v is a solution of the minimization problem. This

completeness of V with respect to · a , and hence with respect to · 1 , is

not satis¬ed in the de¬nition (2.7), as Example 2.8 has shown. Therefore,

an extension of the basic space V , as formulated in (2.20), is necessary.

This space will turn out to be “correct,” since it is complete with respect

to · 1 .

3.1. Variational Equations and Sobolev Spaces 93

In what follows we use the following general assumption:

V is a vector space with scalar product ·, · and the norm ·

1/2

induced by ·, · (for this, v := v, v for v ∈ V issatisf ied) ;

V is complete with respect to · , i.e. a Hilbert space ; (3.1)

a : V — V ’ R is a (not necessarily symmetric) bilinear form ;

b : V ’ R is a linear form .

The following theorem generalizes the above consideration to nonsym-

metric bilinear forms:

Theorem 3.1 (Lax“Milgram) Suppose the following conditions are sat-

is¬ed:

• a is continuous (cf. (2.42)); that is, there exists some constant

M > 0 such that

|a(u, v)| ¤ Mu for all u, v ∈ V ;

v (3.2)

• a is V -elliptic (cf. (2.43)); that is, there exists some constant ± > 0

such that

a(u, u) ≥ ± u for all u ∈ V ;

2

(3.3)

• b is continuous; that is, there exists some constant C > 0 such that

|b(u)| ¤ C u for all u ∈ V . (3.4)

Then the variational equation (2.21), namely,

¬nd u ∈ V such that for all v ∈ V,

¯ a(¯, v) = b(v)

u (3.5)

has one and only one solution.

Here, one cannot avoid the assumptions (3.1) and (3.2)“(3.4) in general.

Proof: See, for example, [26]; for an alternative proof see Exercise 3.1. 2

Now returning to the example above, the assumptions (3.2) and (3.3) are

obviously satis¬ed for · = · a . However, the “preliminary” de¬nition of

the function space V of (2.7) with norm · a de¬ned in (2.19) is insu¬cient,

since (V, · a ) is not complete. Therefore, the space V must be extended.

Indeed, it is not the norm on V that has been chosen incorrectly, since V

is also not complete with respect to another norm · that satis¬es (3.2)

and (3.3). In this case the norms · and · a would be equivalent (cf.

(2.46)), and consequently,

(V, · ⇐’ (V, · ) complete .

a) complete

Now we extend the space V and thereby generalize de¬nition (2.17).

94 3. Finite Element Methods for Linear Elliptic Problems

De¬nition 3.2 Suppose „¦ ‚ Rd is a (bounded) domain.

The Sobolev space H k („¦) is de¬ned by

H k („¦) := v : „¦ ’ R v ∈ L2 („¦) , the weak derivatives ‚ ± v exist

in L2 („¦) and for all multi-indices ± with |±| ¤ k .

A scalar product ·, · and the resulting norm · in H k („¦) are de¬ned

k

k

as follows:

‚ ± v ‚ ± w dx ,

v, w := (3.6)

k

„¦ ± multi’index

|±|¤k