<< . .

. 13
( : 59)



. . >>

. . .
¬ ·
¬ ·
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
$$$
‰
rrr‰r 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

<< . .

. 13
( : 59)



. . >>