¨
4. Ubungsblatt
Iterative L¨
osung großer Gleichungssysteme, UE (106.082), SS 2015
1. a) Vollziehen Sie den Beweis von Lemma 9.1 aus dem Vorlesungsskriptum nach.
b) Exercise 9.1:
Show by means of an induction argument with respect to powers of A :
p(A) r0 = Vm p(Hm )VmT r0 = kr0 k2 Vm p(Hm ) e1
for all p ∈ Pm−1 , where e1 is the first unit vector in IRm .
2. Exercise 9.4:
Assume that for m = 1, 2, . . ., the tridiagonal matrices
α1 β2
β2 α2
β3
..
..
..
Tm =
.
.
.
βm−1 αm−1 βm
βm αm
∈ IRm×m
from (9.12) admit LU-decompositions Tm = Lm Um . Of course, Lm
adopt the notation
1
η1 ω2
λ2 1
η2
.
.
..
..
Tm = Lm Um =
λ
1
m−1
λm 1
and Um are bidiagonal, and we
ω3
..
.
..
.
ηm−1 ωm
ηm
(1)
a) Verify the following recursive formulas for the values λm , ωm , ηm in (1):
ωm = βm ,
λm =
βm
,
ηm−1
ηm = αm − λm ωm
(2)
Conclude that the matrices Lm and Um are recursively obtained from Lm−1 , Um−1 by adding one
row and column, i.e.,
0
0
Um−1
Lm−1
ω
m
,
U
=
(3)
Lm =
m
T
T
0
λm
1
0
η
m
b) Given the factors L, U ∈ IRm×m of the LU-decomposition of a matrix T ∈ IRm×m in the form (3),
L0 0
U0 u
L =
U =
T
T
,
`
1
0
η
with L0 , U 0 ∈ IR(m−1)×(m−1) and `, u ∈ IRm−1 , conclude that L−1 and U −1 can be written as
0 −1 − 1 U 0 −1 u
0 −1
U
L
0
η
U −1 =
L−1 =
,
0T
−1
1
T
0
−` L
1
η
This means that after adding one row and column to the tridiagonal matrix T to obtain the new
tridiagonal matrix T 0 , the triangular inverses L0 −1 , U 0 −1 can be also obtained from L−1 and U −1
by simply adding one row and column.
3. Vollziehen Sie den Abschnitt u
¨ber ‘Block Arnoldi’ aus dem Buch von Y. Saad nach (siehe Anhang),
¨
und erl¨autern Sie das in der Ubung.
Konkret: Weisen Sie nach, dass Algorithm 6.21 die Block-Verallgemeinerung von Algorithm 9.1 aus
dem Skriptum ist, wobei p die Dimension der Bl¨ocke ist. Checken Sie insbesondere die Orthogonalit¨atseigenschaften der Bl¨
ocke Vj und zeigen Sie, dass (6.108) gilt.
Anmerkung: In Algorithmus 6.21 bedeutet ‘V unitary’: V T V = Ip×p . Mit ‘Q-R factorization’ ist die
reduzierte QR-Zerlegung gemeint.
4. Programmieren Sie die Lanczos-Iteration (Algorithm 9.3) in Matlab in der afun - Variante (also ohne
die Matrix A explizit zu verwenden).
5. Die Hessenberg-Matrizen Hm aus der Arnoldi-Iteration werden in verschiedenem Zusammenhang als
Approximationen verwendet. Will man z.B. die Anwendung der Matrix-Exponentialfunktion 1 exp(A)
auf einen Vektor y numerisch approximieren, so kann man dazu
y(t) = exp(A)y ≈ ym := Vm exp(Hm )VmT y = Vm exp(Hm ) βe1 ,
β := kyk2
verwenden (motiviert durch Exercise 9.1).
Sei A die negative definite Matrix, die aus der FD-Approximation von −u00 entsteht (wie im 1DPoisson-Modellproblem). Berechnen Sie exp(A) y f¨
ur irgendeinen Startvektor y mit der oben angegebenen Krylovraum-Methode (unter Verwendung ihres Codes aus Aufgabe 4.), und und protokollieren
Sie den Fehler der Approximation, z.B. in der k · k2 - Norm, f¨
ur m = 1, 2, . . .. Variieren Sie auch die
Dimension von A (z.B. n = 10, 100, 1000).
Zum Vergleich berechnen Sie exp(A) y exakt mithilfe der Spektralzerlegung von A und der diskreten
Sinustransformation [i]dst (vgl. Abschnitt 4.1).
6. In der Praxis werden Krylovraum-Verfahren u
¨ber Residuenkontrolle gesteuert, d.h., im allgemeinen
wird abgebrochen, wenn die Norm des Residuums unter einer vorgegebenen Toleranz liegt. Was ist nun
das ‘Residuum’ der Methode aus Aufgabe 5? Zeigen Sie:
a) Die Approximation ist skalierungsinvariant: F¨
ur tA anstelle von A ergibt sich genau tHm mit dem
zu A geh¨
origen Hm .
b) Damit ist
ym (t) := Vm exp(tHm ) βe1 = Vm um (t)
eine ‘kontinuierliche’ Approximation f¨
ur exp(tA), zumindest rein gedanklich. Die Koeffizientenfunktion um (t) ist L¨
osung des Anfangswertproblems
u0m (t) = Hm um (t),
um (0) = βe1 .
Zeigen Sie: Das Residuum (auch: ‘Defekt’)
0
rm (t) = ym
(t) − A ym (t)
von ym (t) bez¨
uglich der Differentialgleichung x0 = A x ist gegeben durch
rm (t) = −hm+1,m eTm exp(tHm ) β e1 vm+1
|
{z
}
um
und ist daher orthogonal zu span{v1 , . . . , vm }.
Anmerkung: Dieses Residuum ist ein Maß f¨
ur die G¨
ute der Approximation. In der Literatur wird
vorgeschlagen, dieses an mehreren Stellen t auszuwerten und daraus durch grobe Approximation der
Differentialgleichung f¨
ur den Fehler em (t) = ym (t) − y(t) = ym (t) − exp(tA) y,
e0m (t) = A em (t) + rm (t),
em (0) = 0
den Fehler zu sch¨
atzen.
Kleine Zusatzaufgabe: Wenn Sie mit einem ‘glatten’ Vektor y starten – wie sehen die Krylov-Basisvektoren
vi , i = 1, 2, 3, . . . aus? Plotten Sie.
1
bzw. irgendeine Matrixfunktion f (A)
³
+ µ#Q| ` j
|¢V|G#O*
Z+#
$ and an approximation are
; :
] 0
0 M
M
O O
] 0
0 M
$ M
$ O /O 0
Y
An explicit expression for the coefficient readily obtained from (6.99-6.100) by taking
$ ¬
M
$ ? M
? M
6
:
$ M
²
²
²
¬
²
²
²
Since the condition number
of the matrix of eigenvectors is typically not
known and can be very large, results of the nature of the corollary are of limited practical
6
interest. They
can: be useful only when it is known that the matrix is nearly normal, in
which case,
.
O
~E }E ¦ V
¡'£ In many circumstances, it is desirable to work with a block of vectors instead of a single
vector. For example, out-of-core finite-element codes are more efficient when they are
programmed to exploit the presence of a block of the matrix in fast memory, as much as
possible. This can be achieved by using block generalizations of Krylov subspace methods,
for which always operates on a group of vectors instead of a single vector. We begin by
describing a block version of the Arnoldi algorithm.
§
L 3~
% §
)+"! * ,.-'/
}¬
&$%,2* 23
¤¥
1. Choose a unitary matrix
of dimension
2. For c bWYWYZY8 Do: 9
O
4 K K 3.
Compute 4 K
bWYZYWY88
O4 4 K
K
4
4.
Compute K
9 of K : K
5.
Compute the Q-R factorization
6. EndDo
¬ § ¬
¬ § ²
¬
.
K
79
K
7\9
K
The above algorithm is a straightforward block analogue of Algorithm 6.1. By construction, the blocks generated by the algorithm are orthogonal blocks that are also orthog9 identity matrix and use the
onal to each other. In the following we denote by the 9
0
following notation:
¬
¬
¥
¬
6
WYZYWYW &
9 4K : 4 K 9 matrix of the last 4K <
;
for - M
columns of Y
«
O
³
j \ ¦V|G#O
6
+ ® TiM° =±
Then, the following analogue of the relation (6.4) is easily proved:
§
¬
M
79
79 Y
Here, the matrix is no longer Hessenberg, but band-Hessenberg, meaning that it has
subdiagonals instead of only one. Note that the dimension of the subspace in which the
solution is sought is not but Y .
A second version of the algorithm uses a modified block Gram-Schmidt procedure
instead of the simple Gram-Schmidt procedure used above. This leads to a block generalization of Algorithm 6.2, the Modified Gram-Schmidt version of Arnoldi’s method.
L D~
L¬
& $%,2* 23
) ! * ,.-'/
¤¥
1. Choose a unitary matrix
of size
c bZ YWYZY8! Do: 9
2. For O
K
3.
Compute K
4.
For 4LK O! bW4ZY YWY38
K do:
5.
K ! K
4 4K
6.
7.
EndDo
8.
Compute the Q-R decomposition K
9. EndDo
v¬
¬
¬ §
¬
* ,
35
m
-"/
²
¬
K
79
K
79 K
Again, in practice the above algorithm is more viable than its predecessor. Finally, a third
version, developed by A. Ruhe [170] for the symmetric case (block Lanczos), yields a variant that is quite similar to the original Arnoldi algorithm. Assume that the initial block of
orthonormal vectors, ZYWYWYZ is available. The first step of the algorithm is to multiply
9 the resulting
:
by and orthonormalize
vector against ZYWYZYW . The resulting vector
9is defined to be 9
:
. In the second step it is that is multiplied by and orthonormalized
W
:
7
9
against all available 4 ’s. Thus, the algorithm works similarly to Algorithm 6.2 except for
a delay in the vector that is multiplied by at each step.
§
L D~
) ! * ,.-'/
& $%,2* 23
§
)
§
% # "('&36( $5
1. Choose initial orthonormal vectors 4 4
.
9 :
M
W
Z
Y
W
Y
8
Y
2. For Do:
!
O M ;
3.
Set 9
!
4.
Compute O ;
6
5.
For bWYZYWY38
: 0 Do:
!
O
4
6.
4
!0
4 4
7.
0
8.
EndDo
! !
K
9.
Compute and K
.
K
7
9
0
7
9
7 9 0
\
10. EndDo
L¬
v¬
¬ z²
¬
¬ ²
¬ §
¬
¬
¬
Observe that the particular case coincides with the usual Arnoldi process. Also, the
O
dimension of the subspace of approximants, is no longer restricted to being a multiple
© Copyright 2024 ExpyDoc
