Numerische Algorithmen der linearen Algebra - pg.gda.planiem/dyd-zips/Vorlesung-2-beamer.pdf ·...
49
Numerische Algorithmen der linearen Algebra A. Niemunis IBF-Karlsruhe Karlsruhe, 7. November 2018
Transcript of Numerische Algorithmen der linearen Algebra - pg.gda.planiem/dyd-zips/Vorlesung-2-beamer.pdf ·...
Numerische Algorithmen der linearen Algebra
A. NiemunisIBF-Karlsruhe
Karlsruhe, 7. November 2018
Qualitaten der Matrizen (1)
I quadratische, symmetrische M.
I semipositiv- positiv- negativ-definite M.
I diagonal-dominante M., bi-, tri- , diagonale M. DiagonalMatrix[{1,2,3}]
I ahnliche M. A u. B erfullen A = S−1 ·B · SI sparse (= schwach besetzte) / volle Matrix
I dicht gepackte SparseArray[{{1,1} -> 1, {2,2}-> 2, {1,3}-> 4}] // Normal
I dreieckige M., Block-M.
I gut konditionierte M. LinearAlgebra‘MatrixConditionNumber[{{10^6,0},{1.0,10^-6}}]
I orthogonale M. ≈ unitare M. Inverse[ Transpose[A] ] - A; Orthogonalize[..]
I voll, band- M. SparseArray[{Band[{1,1}]-> x, Band[{2,1}]->y}, {5,5}]
I Bildraum ≈ Spalentraum bilden alle y aus y = A · x.Wenn A singular ist, dann y = A · x hat mit y im Bildraum eine mehrdeutige- sonst
keine Losung x.
I (Zeilen)Rang =Anzahl der Lin-unabhang.Zeilen MatrixRank[ m ]
I Kern = bilden alle x-Vektoren mit A · x = 0 NullSpace[{{1,2,3},{4,5,6}}]
I Defekt = Dimension des Kerns
Lineare Gleichungssysteme AN×N · xN×1 = bN×1 (2)
Gauß Algorithmus: Losung beliebiger linearer Gleichungssysteme.Elimination: aquivalentes Gleichungssystem mit dreieckiger Matrix A durchManipulation mit den Zeilen.
A b
z1
z5z6z7
z6z7
Der folgende Schritt eliminiert diefolgenden KomponentenA54 → 0, A64 → 0, A74 → 0
z5 = z5 − z4A54
A44
z6 = z6 − z4A64
A44
z7 = z7 − z4A74
A44
Voraussetzung: A44 6= 0, sonst Pi-votisierung.Achtung: Zeile beinhaltet auch b.
Die Null-Spalten links unten werden dadurch nicht zerstort.
Gauß Elimination (2a)
GaussElimination [aIn_, Id_ ] := Module[{a, d, ud, j, m, n},
a = aIn; {n, m} = Dimensions[a];
If[Id <= 0 || Id >= n, Print["Id = ", Id, " is out of range"]];
d = a[[ Id, Id]];
Do[ ud = a[[ j, Id ]]; a[[j, All ]] = a[[j, All ]]-a[[Id, All ]]*ud/d ;, {j, Id + 1, n}];
a ]
a = Array[ -5 + 10*Random[] & , { 7, 7} ]; a += Transpose[ a] ;
g = Array[0 &, 7]; g[[1]] = MatrixPlot[ a, ColorFunction -> "DarkRainbow"] ; b=a;
Do[b = GaussElimination [b , i ]; g[[i + 1]] = MatrixPlot[ b, ColorFunction -> "DarkRainbow"],
{i, 1, 6}]; g
1 2 3 4 5 6 7
1
2
3
4
5
6
7
1 2 3 4 5 6 7
1
2
3
4
5
6
7
1 2 3 4 5 6 7
1
2
3
4
5
6
7
1 2 3 4 5 6 7
1
2
3
4
5
6
7
1 2 3 4 5 6 7
1
2
3
4
5
6
7
1 2 3 4 5 6 7
1
2
3
4
5
6
7
1 2 3 4 5 6 7
1
2
3
4
5
6
7
1 2 3 4 5 6 7
1
2
3
4
5
6
7
1 2 3 4 5 6 7
1
2
3
4
5
6
7
1 2 3 4 5 6 7
1
2
3
4
5
6
7
1 2 3 4 5 6 7
1
2
3
4
5
6
7
1 2 3 4 5 6 7
1
2
3
4
5
6
7
1 2 3 4 5 6 7
1
2
3
4
5
6
7
1 2 3 4 5 6 7
1
2
3
4
5
6
7
Lineare Gleichungssysteme AN×N · xN×1 = bN×1 (3)
Ruckwartseinsetzen: Die aquivalente dreieckige Form desGleichungssystems lasst sich in der Reihenfolge von unten nach oben losen:
A bx
=
=
.
.
Die Unbekannten: x7 dann x6
dann x5 . . . werden sequentiell be-rechnet
x7 = b7/A77
x6 = (b6 −A67x7)/A66
x5 = (b5 −A57x7 −A56x6)/A55
. . .
Ruckwartseinsetzen (3a)
1 2 3 4 5 6 7 8
1
2
3
4
5
6
7
1 2 3 4 5 6 7 8
1
2
3
4
5
6
7
1 2 3 4 5 6 7 8
1
2
3
4
5
6
7
1 2 3 4 5 6 7 8
1
2
3
4
5
6
7rh
s
rhs'
GaussBacksubstitution [aIn_ ] := Module[{a, b, m, n, ai, x },
{n,m}=Dimensions[aIn]; If[ m !=n+1,Print[ "no RHS or multiple RHS"]];
a = aIn[[All, 1 ;; n]]; b = aIn[[All, n + 1]]; x = Array[0 &, n];
Do[ai= a[[i,All]]; x[[i]]=(b[[i]] - ai. x )/a[[i, i]], {i, n,1,-1 }];
x
]
a = Array[ -5 + 10*Random[] & , { 7, 7} ]; a += Transpose[ a] ;
rhs = Array[-5 + 10*Random[]&,7]; b=Append[a, rhs]// Transpose;
{m,n} = Dimensions[b;] Do[b= GaussElimination[b,i],{i,1,n-1}];
GaussBacksubstitution [b ]
LinearSolve[a, rhs] (* check it *)
Lineare Gleichungssysteme AN×N · xN×1 = bN×1 (4)
Symbolische Notation fur Elimination und Pivotisierung
I Zeile i wird mit einem Faktormultipliziert und zur Zeile jaddiert um Aji = 0 zubekommen. Formellmultiplizieren wir A · Gi mit dersog. Frobeniusmatrix Gi.
I Zeilen (Spalten) i, j vertauschenwir um stabiles Pivotelement 6= 0zu wahlen. Formell erfolgt diesdurch Pra-(Nach)-multiplikationmit Permutationsmatrix Pi−j
z.B. P1−4 vertauscht z1 und z4 →
G1 =
1 0 0 0
−a21/a11 1 0 0
−a31/a11 0 1 0
−a41/a11 0 0 1
P1−4 =
0 0 0 1
0 1 0 0
0 0 1 0
1 0 0 0
Ruckgangig mit: (Gi)
−1 = 2 1− Gi und P−1i−j = Pi−j
Pivotisierung (4a)
Pivotisierung = geschickter Tausch der Zeilen in A und in b um dieNull-Division im Gauß-Algorithmus zu meiden.
Ohne Pivotisierung konnte der Gauß Algorithmus fur A · x = b mit derMatrix
A =
0 1 2 3
4 5 6 7
. . . . . . . . . . . . 1
. . . . . . . . . . . . 1
nicht mal starten, da schon bei der ersten Elimination
z2 = z2 − z1A21
A11
eine Null-Division auftaucht.
Nach der Losung muss der Tausch ruckgangig gemacht werden.
Lineare Gleichungssysteme: LU- = LR-Zerlegung (5)
AN×N unsymm.: A = L · U (oder A = P · L · U mit Pivotisierung)Die LU-Zerlegung der Matrix A mit dem Crout’schen Algorithmus
A = L · U =
∗ 0 0 0
∗ ∗ 0 0
∗ ∗ ∗ 0
∗ ∗ ∗ ∗
·
∗ ∗ ∗ ∗
0 ∗ ∗ ∗
0 0 ∗ ∗
0 0 0 ∗
L · U · x = b wird via Vorwartseinsetzen L · c = b gelost und mit demZwischenergebnis c folgt die Ruckwartseinsetzen U · x = c.
Beim Gauß- wird aus der Matrix A die Matrix U berechnet. Die Komponenten von L errechnensich automatisch als die dabei verwendeten Zeilen-Multiplikatoren (wie in der FrobeniusmatrixG), z.B. L21 = A21/A11, L31 = A31/A11 usw.
Lineare Gleichungssysteme: LU- = LR-Zerlegung (5)
AN×N unsymm.: A = L · U (oder A = P · L · U mit Pivotisierung)Die LU-Zerlegung der Matrix A mit dem Crout’schen Algorithmus
A = L · U =
1 0 0 0
∗ 1 0 0
∗ ∗ 1 0
∗ ∗ ∗ 1
·
∗ ∗ ∗ ∗
0 ∗ ∗ ∗
0 0 ∗ ∗
0 0 0 ∗
L · U · x = b wird via Vorwartseinsetzen L · c = b gelost und mit demZwischenergebnis c folgt die Ruckwartseinsetzen U · x = c.
Bei Gauß-Elim. wird aus der Matrix A die Matrix U berechnet. Die Komponenten von Lerrechnen sich automatisch als die dabei verwendeten Zeilen-Multiplikatoren (wie in derFrobeniusmatrix G), z.B. L21 = A21/A11, L31 = A31/A11 usw.
LU- = LR-Zerlegung mit M-ma, Beispiel (5a)
{LU,p,c}= LUDecomposition[{{1,2,3}, {1,1,1}, {3,3,1}}]
{{{1, 2, 3}, {1, -1, -2}, {3, 3, -2}}, {1, 2, 3}, 1}
LU =
[1 2 31 −1 −23 3 −2
]besteht aus Uberlagerung von 2 Matrizen: LU = L+ U − I3x3, d.h.
L =
[0 0 01 0 03 3 0
]+
[1 0 00 1 00 0 1
]und U =
[1 2 30 −1 −20 0 −2
]p zeigt die als Pivot verwendeten Zeilenc zeigt die Konditionierungszahl der Matrix
Lineare Gleichungssysteme Determinante (6)Eigenschaften und Berechnung von detA
I 2 Zeilen vertauschen ⇒ −1 ∗ det, z.B. Pivotisierung
I mehrfaches Addieren einer Zeile /Spalte ⇒ 1 ∗ det
I Multiplizieren einer Zeile/Spalte mal λ⇒ λ ∗ det.Sonderfall: det(λA) = λN det(A)
I Transposition ⇒ 1 ∗ det
I det(A · B) = det(A) det(B)
I Dreieckige Block-Matrix⇒ det(A) = det(A11) det(A22) . . . .
xe
e
x
1
2
1
2
det[ ]
det[ ]u
u
v
vw
- u - - v -- w -
- u - - v -
A
A
A
A
A
11
22
33
44
55
0
Sonderfall: detU = a11a22 . . . aNN
a
aa
a
a
a
11
22
33
44
55
66
det(L) det(U) = detU , mit det(L) = 1
Zu Block-Matrizen siehe auch http://en.wikipedia.org/wiki/Schur_complement
Nutzliche Operationen mit Block Matrizen (6a)
Multiplikation von Blockmatrizen mit quadratischen diagonalen Blockengleicher Dimension erfolgt analog wie bei skalaren Komponenten:[
A1 B1C1 D1
]·[A2 B2C2 D2
]=[A1 ·A2 +B1 · C2 A1 ·B2 +B1 ·D2C1 ·A2 +D1 · C2 C1 ·B2 +D1 ·D2
](1)
Inversion der blockdiagonalen Matrix:[A
D
]−1
=
[A−1
D−1
]und det
[A
D
]= detAdetD (2)
Faktorisierung beim nicht-singularen Block[A BT1B2 C
]=[
I 0B2A
−1 I
]·[A 00 C −B2A
](3)
Inversion einer Matrix A→ A−1 (7)
Der Gauß-Jordan Algorithmus: LU-Zerlegung und b→ B = {e1, e2, . . . }.ei sind Achsenvektoren, z.B. e3 = {0, 0, 1, 0}T[ ∗ 0 0 0
∗ ∗ 0 0∗ ∗ ∗ 0∗ ∗ ∗ ∗
]·
[ ∗ ∗ ∗ ∗0 ∗ ∗ ∗0 0 ∗ ∗0 0 0 ∗
]·
[x y z vx y z vx y z vx y z v
]=
1 0 0 00 1 0 00 0 1 00 0 0 1
Die Losungen x,y, z,v der Gleichungen A · x = e1, . . . , A · v = e4
konnen als Spalten der inversen Matrix A−1 gesehen werden
A · A−1 = 1
mit der Identitatsmatrix 1 (6= Einheitsmatrix)
Zeitaufwand bei Inversion von An×n (8)
Gauß VerfahrenGauß Elimination: in etwa n3/3 FLOP s = floating point operation
Gauß Ruckwarts-Substitution: in etwa n2/2 FLOP s
Alles wird n− fach ausgefuhrt, d.h. n4/3 + n3/2 /
LU-VerfahrenL-U Elimination: n3/3 FLOP sL-U Ruckwarts-Substitution: n2/2 FLOP sL-U Vorwarts-Substitution: n2/2 FLOP sNur Substitutionen werden n− fach ausgefuhrt, d.h. : n3/3 + nn2/2 = 4
3n3
FLOP s ,3-fach feineres 2D FE-Netz 1 min → 12h
Lineare Gleichungssysteme (9)
Cholesky-Verfahren (nur fur positive symm. Matrix A)
Zerlegung A = L · LT = (L · D1/2) · (D1/2 · LT ) = L · D · LT moglich fur jedeN ×N Matrix, die positiv definit ist∀x gilt x · A · x > 0.
L =
1 0 0 0∗ 1 0 0∗ ∗ 1 0∗ ∗ ∗ 1
Die Pivotiseirung ist fur positive symmetrische Matrizen nicht notig
Band-Diagonale Gleichungssysteme (10)
Tridiagonale Systeme sind besonders einfach.Bei diagonaler Dominanz erubrigt sich die Pivotisierung.
c6
b6a6 c6
b6a6
Diagonale Dominanz:|bi| > |ai|+ |ci|
Heute ist eine kurze CPU-Zeit wichtiger als Speicherbedarf.Ist Speicherbedarf entscheidend dann frontale Gleichungsloser.
Tridiagonalmatrix fur axial belasteten Pfahl (11)
kp
ksFL
1 2
1 2 1 2uu
1 2
RR11
2 3 4 5 6 7 8
uα = Verschiebung (nach rechts uα > 0) vom Knoten α ;
Reα = Kraft auf Knoten α vom Stabelement e;
Fα = Externe Kraft auf Knoten α.
ks = ksπDL/2 = Feder aus der Bettung uber die Halfte derMantel-Flache des Segments
kp = kpπD2/2 = Feder aus der Bettung unter dem Pfahlfuß
Beachte lokale (innerhalb Segments) und globale Nummerierung
Interne Krafte auf die Knoten � (Stab + Bettung) (12)
jedes Element e hat Knoten 1,2 (lokale Nummerierung)
1 2ε > 0u
u1
1
2
2
u1 u 2
positive Kräfte u. Verschiebungen nach rechts
R
1
1R
R
2
2Rs s
e e
ex Aus dem Stab e (Zug > 0):
R1e = AEε =
EA
L(u2−u1)
R2e = −AEε = −R1
e
Aus der Bettung:
R1s = −ksπD 1
2L u1
R2s = −ksπD 1
2L u2
Zusammen aus dem Element Nr e und aus der Bettung entlang e:[R1R2
]= −
[k + ks −k−k k + ks
]e·[u1u2
]oder R = −Ke · ue
mit Abkurzungen k = EA/L und ks = ksπD12L
Tridiagonalmatrix fur axial belasteten Pfahl (13)
kp
ksFL
1 2
1 2 1 2uu
1 2
RR11
2 3 4 5 6 7 8
Fur Gleichgewicht am Knoten Nr i (globale Nummer) addiert man die Kraft R2
aus dem Element e = i− 1 (links vom Knoten i) und R1 aus dem Element e = i(rechts vom Knoten i) und evtl. die externe Kraft Fi (bei uns nur F1 6= 0), d.h.
∀i : R2e=i−1 +R1
e=i + Fi = 0 mit[R1R2
]= −
[k + ks −k−k k + ks
]e·[u1e
u2e
]Nach diesem Muster:
��ku0 − (k + ks)u1 +ku2 + F1 =0 fur i = 1
ku1 − 2(k + ks)u2 +ku3 +��F2 =0 fur i = 2
ku2 − 2(k + ks)u3 +ku4 +��F3 =0 fur i = 3
. . .
ku6 − 2(k + ks)u7 +ku8 +��F7 =0 fur i = 7
ku7 − [(k + k∗s ) + kp]u8 +��ku9 +��F8 =0 fur i = 8
Das System schreibt man in Matrixform mit Kraften Fi auf der rechten Seite.
Tridiagonalmatrix (diagonal dominant) (14)
a −k 0 0 0 0 0 0u
−k b −k 0 0 0 0 0u
0 −k b −k 0 0 0 0u
0 0 −k b −k 0 0 0u
0 0 0 −k b −k 0 0u
0 0 0 0 −k b −k 0u
0 0 0 0 0 −k b −ku0 0 0 0 0 0 −k cu
·
uu1uu2uu3uu4uu5uu6uu7uu8
=
F10000000
mita = k + ks,b = 2(k + ks)c = k+ ks + kp
Die Struktur ist typisch fur die FEund resultiert aus der Zusammenset-zung der Beitrage zur Steifigkeit ausden einzelnen Elementen.Evtl. Beitrage aus der Streckenlast werden
auch zusammengesetzt zu F.
Ke
e=7
e=1
e=2...6
Kglob
-k
-k
-k
-k
-k
-k
-k
-k
-k
-k
-k
-k
-k
-k
-k
-k
-k
-k
-k
-k
k
-k
-k
k
=
F
b
b
b
b
b
b
cc
a
a
a
a
a
a
Nichtlinearitat (14a)
τ
γ
τ
γ
τ
γ
Nichtlinearitat bei einer starkenaxialen Belastung eines Pfahls(nach Randolph und Wroth 1978).
Eine inkrementelle NL Losung
ist erforderlich. /
A · x = b→ iterative Verbesserung der Losung (15)
Eine schlecht konditionierte Matrix A mit der Zerlegung A = L · U fuhrt zur(ungenauen) numerischen Losung x + δx. Sie erfullt nur naherungsweise dieGleichung
A · (x + δx) ≈ b (4)
d.h., es ergibt sich eine Ungenauigkeit, sog. Residuum∗:
A · (x + δx)− b = δb 6= 0 (5)
Den Fehler δx findet man schnell aus
δx = (L · U)−1 · δb
Die Korrektur −δx wird zur Losung x + δx addiert
x = (x + δx)− δx (6)
∗Streng genommen, als Residuum bezeichnet man −δb in der Literatur
Eigenwerte von unsymm. An×n (16)Eine n× n Matrix A hat n Eigenwerte. Sie ist diagonalisierbar wennsie auch n linke (oder rechte) lin. unabhangige Eigenvektoren hat:
l(i) ·A = λ(i)l(i) bzw. A · r(i) = λ(i)r(i) mit i = 1, 2, . . . n (7)
Das Spektrum Λ = diag(λ(1), λ(2), . . . λ(n)) und linke l(i) oder rechte r(i) Eigenvektoren
konnen komplex sein. Alle l(i) und r(i) sind NORMIERT.
Die Spaltenmatrix R ={r(1), r(2), . . .
}Tist invertierbar, wenn alle
Spalten r(i) lin. unabhangig sind. Aus A ·R = R ·Λ folgt dann dieSpektralzerlegung A = R ·Λ ·R−1.
Die Zeilenmatrix L ={
l(1), l(2), . . .}
ist invertierbar, wenn alle Zeilen
l(i) lin. unabhangig sind. Aus L ·A = Λ · L folgt dann dieSpektralzerlegung
A = L−1 ·Λ · L = R ·Λ ·R−1 (8)
Das Spektrum Λ ist fur die linke und rechte Eienvektoren identisch. Auch die linke- und die
rechte inverse Matrizen sind gleich A ·A−1 = A−1 ·A = I. Fur a 6= b gilt r(a) · l(b) = 0 aber
i.Allg. r(a) · r(b) 6= 0 und l(a) · l(b) 6= 0 und r(a) · l(a) 6= 1. Mit anderen Worten, R und L sind
keine orthogonale Matrizen und R 6= L−1. Eine nutzliche Zerlegung der
nicht-diagonalisierbaren Matrizen erfolgt mit SDV.
Eigenwerte von unsymm. An×n (17)
Die Spektralzerlegung
A =[R ·Λ ·R−1
]=[L−1 ·Λ · L
]=
n∑i
λ(i)r(i)l(i)
r(i) · l(i)
hat die Matrixform
A =
[↑ ↑ ↑r(1) . . . r(n)
↓ ↓ ↓
]·
λ(1) . . .λ(n)
· [ ↑ ↑ ↑r(1) . . . r(n)
↓ ↓ ↓
]−1
=
=
← l(1) →
←... →
← l(n) →
−1
·
λ(1) . . .λ(n)
·← l(1) →
←... →
← l(n) →
Bei normierten Eigenvektoren ‖r(i)‖ = ‖l(i)‖ = 1 gilt R−1 6= L aberbei einer alternativen Skalierung der linken Eigenvektoren
l(i)
= l(i)/(r(i) · l(i)) gilt L = R−1 auf Kosten von ‖l(i)‖ 6= 1.
SVD Zerlegung einer unsymm. Matrix An×n (18)Singular-Value-Decomposition (=SVD)
A = U ·D ·VT mit UT ·U = 1 und VT ·V = 1 singulare Werte = D 6= Λ
Im Gegensatz zur spektralen Zerlegung unterschiedliche Rotationen VT und U (und
Skalierung mit D statt Λ )
Die diagonale Matrix D hat positive Komponenten Di ≥ 0,∈ R
Ein SVD Solver:[U ·D ·VT
]· x = b oder A · x = b
D ·VT · x = UT · bwenn D invertierbar x = (D ·VT )−1UT · b
x =[V ·D−1 ·UT
]· b oder x = A−1 · b
AA = {{0.0, 1}, {-5, 0}};
Lamb = Eigenvalues[AA]; (*complex*)
{Lamb, {v1,v2}} = Eigensystem[AA] (* complex*)
{UU, DD, VV} = SingularValueDecomposition[AA]
(* = {{{0,-1}, {1, 0}}, {{5, 0}, {0, 1}}, {{-1, 0}, {0, -1}}} *)
Direkte lineare Gleichungsloser (19)
Direkte Gleichungsloser:Gauß, LU-Zerlegung, L.D.LT -Zerlegung (symm.), Cholesky =L.LT -Zerlegung (symm. + positiv)
I bei unbegrenzter Prazision ware die Losung exakt,
I numerische Rundungsfehler konnen bei schlecht konditioniertenMatrizen erheblich sein,
I schwache Effizienz bei großen Matrizen
I Algorithmus ist stabil
I kaum Zeitersparnis bei schwach besetzten Matrizen
I unberechenbare Speicherzugriffe (schwierig parallelisierbar)
Iterative lineare Gleichungsloser (20)
Iterative Gleichungsloser:Stationar: Gauß-Seidl, Jacobi, Uber-RelaxationNicht stationar : (Bi-)Konjugierte Gradienten (stabilisiert / quadriert ),(Quasi) Min. Residuum (verallgemeinert, BFGS), Chebyschev IterationBroyden-Fletcher-Goldfarb-Shanno
Multigrid: sehr schnell (algebraisch = ohne Netz)
I Konvergenzkriterium definieren
I gute Effizienz aber instabil
I i.d.R. angepaßt an spezielle Matrizen
I leichter parallelisierbar
Gram-Schmidt Orthonormalisierung (21)
Aufgabe:N linear unabhangige Vekto-ren v1,v2, . . . ,vN in Rn mit(N ≤ n) spannen einen N -dimensionalen Raum auf, sog.span(v1,v2, . . . ,vN ).
v1
u1u2v2
v3u3
v2v1
v3
Span = lineare Hulle. Dazu gehoren alle lineare Kombinationen∑
i αivi.
Aus v1,v2, . . . ,vN finden wir ~u1, ~u2, . . . , ~uN die den gleichen Raumaufspannen aber senkrecht zueinander stehen.Gram-Schmidt Orthonormalisierung wird oft fur Eigenvektoren aus einem mehrfachen
Eigenwert benutzt.
Projektion vom Vektor v auf dieWirkungslinie des Vektors u:
proju(v) = (v · ~u) ~u =(v · u)u
u · u
uu v
v
u
v
(u v)u
u
u
v
proj (v) =
Gram-Schmidt Orthonormalisierung (22)
Orthogonalisierung = Berechnung der Reste nach den Projektionen.
u
v
u
u
v v
proj (v)
v- uproj (v)
u1 = v1
u2 = v2 − proju1(v2)
u3 = v3 − proju1(v3)− proju2
(v3)
u4 = v4 − proju1(v4)− proju2
(v4)− proju3(v4)
. . .
Z.B. Vektor v3 vermindert um seine Projektionen auf u1 und u2
ergibt den Rest u3, der senkrecht zu u1 und u2 stehen muss, d.h.
u3 ⊥ u1 und u3 ⊥ u2
QR-Zerlegung einer reellen Matrix AN×N (23)
QR-Zerlegung wird bei Eigenwertproblemen angewandt
A = Q ·R mit QT ·Q = 1
Q ist orthogonal und R ist die rechte obere Dreiecksmatrix.
Verschiedene Algorithmen fur die QR Zerlegung:
I Householder Transformationen
I Givens Rotationen
I Gram-Schmidt’sche Orthogonalisierung
Vorbereitung 1: Gram-Schmidt Verfahren (23a)
Gram-Schmidt Orthogonalisierung der Spalten ai der reellen Matrix AN×N
A = [ a1 | a2 | . . . | aN ]
[ a1|a2| . . . |aN ] → Gram-Schmidt → [ ~u1|~u2| . . . |~uN ] = Q mit
u1 = ~u1‖u1‖ = a1 Euklidische Norm ‖x‖ =
√x · x
u2 = ~u2‖u2‖ = a2 − ~u1~u1 · a2
u3 = ~u3‖u3‖ = a3 − ~u1~u1 · a3 − ~u2~u2 · a3
. . .
Alle ai und ~ui mit i = 1, 2, . . . Nsind nun bekannt unda1 = ~u1‖u1‖a2 = ~u2‖u2‖+ ~u1~u1 · a2
a3 = ~u3‖u3‖+ ~u1~u1 · a3 + ~u2~u2 · a3
. . .
Jonglieren mit Zeilen und Spalten (24)
[ a b c ] ·[ ← x →← y →← z →
]=
[ ↑ ↑ ↑x y z↓ ↓ ↓
]·[ abc
]=
[↑
ax + by + cz↓
]
[ abc
]·[ ← x →← y →← z →
]=
[ ← ax →← by →← cz →
]und
[ ↑ ↑ ↑x y z↓ ↓ ↓
]·[ a
bc
]=
[↑ ↑ ↑ax by cz↓ ↓ ↓
]und
[ ↑ ↑ ↑x y z↓ ↓ ↓
]·[a dbc
]=
[↑ ↑ ↑ax by + dx cz↓ ↓ ↓
]und
[ ↑ ↑ ↑x y z↓ ↓ ↓
]·[a d eb fc
]=
[↑ ↑ ↑ax by + dx cz + ex + fy↓ ↓ ↓
]
QR-Zerlegung (25)Aus dem Gleichungssystema1 = ~u1‖u1‖a2 = ~u2‖u2‖+ ~u1(~u1 · a2)
a3 = ~u3‖u3‖+ ~u1(~u1 · a3) + ~u2(~u2 · a3)
. . .
a1 = a~u1
a2 = b~u2 + d~u1
a3 = c~u3 + e~u1 + f~u2
. . .verglichen mit dem Muster
↑ ↑ ↑
x y z
↓ ↓ ↓
·a d e
b f
c
=
↑ ↑ ↑
ax by + dx cz + ex + fy
↓ ↓ ↓
erkennen wir A = Q ·R
↑ ↑ ↑ ↑
a1 a2 a3 . . .
↓ ↓ ↓ ↓
=
↑ ↑ ↑ ↑
~u1 ~u2 ~u3 . . .
↓ ↓ ↓ ↓
·
‖u1‖ (~u1 · a2) (~u1 · a3) . . .
0 ‖u2‖ (~u2 · a3)
0 0 ‖u3‖...
. . .
Hausholder Transformation (26)
0
u a
a
Spiegelung durch Ebene ⊥ ~u, die 0 beinhaltet:
a = H · a = a− 2~u(~u · a)
mit H = 1− 2~u~u oder Hij = δij − 2~ui~uj
Eigenschaften der Hausholder’schen Matrix:
H = HT = H−1
Aufgabe: Vorgegeben sei Vektor a. Finde H bzw. ~u, das ein gegebenes azur Form a = ρe1 = {ρ, 0, 0, . . . , 0} spiegelt.
‖H · a‖ = ‖a‖ → ρ = ±‖a‖, da H orthog.
a = H · a = a− 2~u~u · a = a− β~u → ~u = (a− a)→
daher: ~u = (ρe1 − a)→
Hausholder Transformation: Anwendung 1 (27a)
QR -Zerlegung einer unsymm. Matrix A[ ∗ ∗ ∗ ∗∗ ∗ ∗ ∗∗ ∗ ∗ ∗∗ ∗ ∗ ∗
]H1·A−→
[ ∗ ∗ ∗ ∗0 ∗ ∗ ∗0 ∗ ∗ ∗0 ∗ ∗ ∗
]H2·H1·A−→
[ ∗ ∗ ∗ ∗0 ∗ ∗ ∗0 0 ∗ ∗0 0 ∗ ∗
]H3·H2·H1·A−→
[ ∗ ∗ ∗ ∗0 ∗ ∗ ∗0 0 ∗ ∗0 0 0 ∗
]
Fur die 1. Spalte :H1 = 14×4 − 2~u~u mit ~u = ({ρ, 0, 0, 0} − {∗, ∗, ∗, ∗})→ und ρ = ‖{∗, ∗, ∗, ∗}‖
Fur die 2. Spalte:H2 = 14×4 − 2~u~u mit ~u = ({0, ρ, 0, 0} − {0, ∗, ∗, ∗})→ und ρ = ‖{0, ∗, ∗, ∗}‖
H2 =
11
11
− 2 ·[
0∗∗∗
]⊗ [ 0 ∗ ∗ ∗ ] =
1 0 0 00 ∗ ∗ ∗0 ∗ ∗ ∗0 ∗ ∗ ∗
(9)
Fur die 3. Spalte:H3 = 14×4 − 2~u~u mit ~u = ({0, 0, ρ, 0} − {0, 0, ∗, ∗})→ und ρ = ‖{0, 0, ∗, ∗}‖
H3 ·H2 ·H1 ·A = R oder A = Q ·R mit Q = (H3 ·H2 ·H1)−1 = H1 ·H2 ·H3
Hausholder Transformation: Anwendung 1 (27b)
getHausholderH[vv_, ee_ ] := Module[{uu, r, n },
n = Dimensions[vv][[1]]; r = Sqrt[vv.vv] // N; uu = Normalize[ (r *ee - vv) ] ;
IdentityMatrix[n] - 2 Outer[Times, uu, uu] ]
HausholderQR[A_] := Module[{B, n, v, H, Hall, e},
B = A; n = Dimensions[B][[2]]; Hall = IdentityMatrix[n];
Do[ v = B[[All, i]]; If[i == 2, v[[1]] = 0 ]; If[i >= 3, v[[ 1 ;; i - 1 ]] = 0 ] ;
e = Array[0 &, n]; e[[i]] = 1;
H = getHausholderH[ v , e ] ; B = H.B; Hall = H.Hall ;
, {i, 1, n - 1}];
{Transpose[Hall], B} ]
Wir uberprufen die oberen Module mit
SP = Array[Random[]&, {100,100}]; {Q, R} = HausholderQR[SP] ; SP-Q.R (* == 0... *)
Hausholder Transformation: Anwendung 2 (27)
Transformation zur bidiagonalen Form ist der erste Schritt derSVD-Zerlegung
A1 =
[ ∗ ∗ ∗ ∗∗ ∗ ∗ ∗∗ ∗ ∗ ∗∗ ∗ ∗ ∗
]A1=H1·A1−→
[ ∗ ∗ ∗ ∗0 ∗ ∗ ∗0 ∗ ∗ ∗0 ∗ ∗ ∗
]A2=A1·H1−→
∗ ∗ 0 00 ∗ ∗ ∗0 ∗ ∗ ∗0 ∗ ∗ ∗
A2=H2·A2−→
∗ ∗ 0 00 ∗ ∗ ∗0 0 ∗ ∗0 0 ∗ ∗
A3=A2·H2−→
∗ ∗ 0 00 ∗ ∗ 00 0 ∗ ∗0 0 ∗ ∗
A3=H2·A3−→
∗ ∗ 0 00 ∗ ∗ 00 0 ∗ ∗0 0 0 ∗
Fur die 1. Spalte :H1 = 14×4 − 2~u~u mit ~u = ({ρ, 0, 0, 0} − {∗, ∗, ∗, ∗})→ und ρ = ‖{∗, ∗, ∗, ∗}‖
Fur die 1. Zeile ohne a11
H1 = 14×4 − 2~u~u mit ~u = ({0, ρ, 0, 0} − {0, ∗, ∗, ∗})→ und ρ = ‖{0, ∗, ∗, ∗}‖Damit bleibt die erste Spalte nach A · H1 intakt.
. . .
Givens Rotation der symm. Matrix A (28)
Givens Matrix Gpq ist orthogonal und hat die Form (hier G36 fur A8×8)
Gpq =
1
1c s
11
−s c1
1
c = cos θ, s = sin θ
GR ist durch GpqT ·A·Gpq defi-
niert. Diese Multiplikation ent-spricht einer aktiven Rotationum Winkel θ in der ”Ebene”xp, xq.
Vorteile:
I Schneller Multiplikationsalgorithmus bei GR
I Der Winkel θ kann so gewahlt werden, dass die Komponenten A∗pq = 0und A∗qp = 0 nach der GR (d.h. A∗ = Gpq
T ·A ·Gpq) verschwinden.
I GR ist gut parallelisierbar
I GR hat kleine Rundungsfehler
Givens Rotation (29)Eine schnelle Multiplikation fur GT
36 ·A ·G36. Nur die Spalten ai3, ai6 unddie Zeilen a3i, a6i von A werden betroffen.
GT36 · A · G36
1
1c −s
11
s c1
1
·
1 ∗ ∗1 ∗ ∗
∗ ∗ a33 ∗ ∗ a36 ∗ ∗1 ∗ ∗1 ∗ ∗
∗ ∗ a36 ∗ ∗ a66 ∗ ∗1 ∗ ∗1 ∗ ∗
·
1
1c s
11
−s c1
1
=
1 • •1 • •
• • • • • 0 • •1 • •1 • •
• • 0 • • • • •1 • •1 • •
Diagonalisierung mit Givens Rotationen (30)
Inkrementelle Diagonalisierung:
1. Finde die betragsmaßig großte außerdiagonale Komponente, Apq
2. Bei max |Apq| < 10−8 ist die Diagonalisierung abgeschlossen.
3. Finde GR fur Apq = Aqp → 0 und rotiere A := GTpq ·A ·Gpq
4. Aktualisiere die globale Rotationsmatrix G := G ·Gpq
Beispiel: Symm. 4× 4 Matrix nach mehreren G-Rotationen.
GR GR
Daraus errechnet sich die spektrale Zerlegung
A = . . .G24 ·G13 · λ ·GT13 ·GT
24 · · · = G · λ ·GT
mit G = G24G13 . . .
Diagonalisierung mit Givens Rotationen (30a)
givensMax [aIn_ ] := Module[{a, p, q, c, s, t, theta, app, apq, aqp, aqq , auxp, auxq},
a = aIn; a = (a - DiagonalMatrix[ Diagonal[a] ]) // Abs ;
{p, q} = Position[ a , Max[a] ] [[1]] // Sort ;
{app, apq , aqp, aqq } = { aIn[[p, p]], aIn[[p, q]], aIn[[q, p]] , aIn[[q, q]]};
If[Abs[apq - aqp] > 10^-15 , Print["Unsymmetric input matrix"]];
If[Abs[apq] < 10^-15, Print["done"]];
theta = (aqq - app) / (2*apq); (* NR 11.1 *)
t = Sign[theta] /(Abs[theta] + Sqrt[theta*theta + 1 ]);
c = 1 / Sqrt[t*t + 1]; s = t*c;
auxp = c*aIn[[All, p ]] - s* aIn[[ All, q ]]; (* new row / column p *)
auxq = c*aIn[[All, q ]] + s* aIn[[ All, p ]]; (* new row / column q *)
a = aIn;
a[[p, All]] = auxp ; a[[All, p]] = auxp ;
a[[q, All]] = auxq; a[[All, q]] = auxq;
a[[p, q]] = 0; a[[q, p]] = 0; (* special cases *)
a[[p, p]] = c*c *app + s*s *aqq - 2* s*c *apq ;
a[[q, q]] = s*s *app + c*c *aqq + 2* s*c *apq ; a
]
a = Array[ -5 + 10*Random[] & , { 10, 10} ]; a = a + Transpose[a];
ev = Eigenvalues[a] // Sort ; (* for check *)
Do[ a = givensMax[a] , {i, 1, 50} ];
MatrixPlot[ a, ColorFunction -> "DarkRainbow"]
Diagonal[a] //Sort (* compare with ev *)
1 2 3 4 5 6 7 8 9 10
1
2
3
4
5
6
7
8
9
10
1 2 3 4 5 6 7 8 9 10
1
2
3
4
5
6
7
8
9
10
1 2 3 4 5 6 7 8 9 10
1
2
3
4
5
6
7
8
9
10
1 2 3 4 5 6 7 8 9 10
1
2
3
4
5
6
7
8
9
10
1 2 3 4 5 6 7 8 9 10
1
2
3
4
5
6
7
8
9
10
1 2 3 4 5 6 7 8 9 10
1
2
3
4
5
6
7
8
9
10
1 2 3 4 5 6 7 8 9 10
1
2
3
4
5
6
7
8
9
10
1 2 3 4 5 6 7 8 9 10
1
2
3
4
5
6
7
8
9
10
LAPACK (31)
LAPACK = Linear Algebra PACKkage in Fortran 77.LAPACK = EISPACK+LINPACK+. . .LAPACK95 = Schnittstelle zu LAPACK in Fortran 95
1. Linear Equations (ES)
2. Linear Least Squares (LLS) Problems
3. Generalized Linear Least Squares (LSE and GLM) Problems
4. Standard Eigenvalue and Singular Value Problems
I Symmetric Eigenproblems (SEP)I Nonsymmetric Eigenproblems (NEP)I Singular Value Decomposition (SVD)
5. Generalized Eigenvalue and Singular Value Problems
I Generalized Symmetric Definite Eigenproblems (GSEP)I Generalized Nonsymmetric Eigenproblems (GNEP)I Generalized Singular Value Decomposition (GSVD)
Hallo Fortran 95 (32)
Fortran 95-Kompiler gfortran ist frei.
Aus
(((((((
(((((((
((((
http://gcc.gnu.org/wiki/GFortranBinaries herunterladen.
Aus meiner Homepage herunterladen und installieren.Erstellen (mit einem Texteditor z.B. mit Notepad++) die Quelltext-Dateiz.B. C:\tmp\nic.f95 mit dem Inhalt
program nic
write(*,*) ’multipliziere 15*16 = ’, 15*16
end program nic
Aus dem DOS-Fenster kompilieren C:\tmp> gfortran nic.f90 -o nic.exe
Aus dem DOS-Fenster ausfuhren C:\tmp> nic.exe
Uber Fortran 90, 95, 2003, 2008 im Internet lesen z.B.http://docs.hp.com/en/B3908-90002/index.html, evtl. GNU Debugger unterhttp://www.gnu.org/software/ddd/
Hallo Lapack95 (33)
C:\pfad\lib einlegen.Dorthin die folgendenDateien1 speichern
blas.a
lapack.a
lapack95.a
tmglib.a
lapack95_modules\f77_lapack.mod
lapack95_modules\f95_lapack.mod
lapack95_modules\la_auxmod.mod
lapack95_modules\la_precision.mod
Lapack und Lapack95 sinddamit verfugbar.
Erstellen C:\tmp\nicL.f90 = unsere Programm
PROGRAM nicL
USE LA_PRECISION, ONLY: WP => DP
USE F95_LAPACK , ONLY: LA_GESV
real( WP ) :: A(2,2), B(2);
A = 1.0_WP*reshape((/ 1,2,4,6/),(/2,2/))
B=1.0_WP*(/1,1/) ;
call LA_GESV( A, B ) ! general equation solver
write(*,*) B
END PROGRAM nicL
Kompilieren mitC:\tmp> gfortran nicL.f90 c:\pfad\lib\lapack95.a c:\pfad\lib\lapack.a
-Ic:\pfad\lib\lapack95_modules -o nicL.exe
Ausfuhren mit C:\tmp> nicL.exe Die Losung x derGleichung A · x = b wird in die Variable B (ur-sprunglich = b) geschrieben und auf dem Bildschirmausgegeben.
Dokumentation der Lapack95-Funktionen liegt unterhttp://www.netlib.org/lapack95/lug95/
(34) Warum ”95”?
Viel FE-software wird in Fortran geschrieben. Fortran95 enthaltFortran (77) aber bietet Erweiterungen, u.a.:
I modules mit public und private
I definierbare Operatoren z.B. if(a .approx. b ) then ...
oder T33 .xx. dev(T33)
I Funktionen, die Matrizen erwidern konnen
I optionale Paremeter
I Polymorphismus z.B. norm(T33)*norm(T3333)
I dynamische Objekte read*,n; allocate( A(n,n))
I eingebaute Funktionen fur Manipulation mit Matrizenz.B. T3 = T44(2:, 2)
I . . .
Fortran2003 bietet zusatzlich Konzepte der OOP (object orientedprogramming), z.B. Klassen, Erbschaft usw.Gfortran wird immer noch zu Fortran2003 entwickelt.Software Sammlung http://www.geotechnicaldirectory.com/cat/Software.html
Numerical Recipes (35)
NUMERICAL RECIPES - The Art of Scientific ComputingWilliam H. PressSaul A. TeukolskyWilliam T. VetterlingBrian P. Flannery
PDF: http://www.nrbook.com/a/bookcpdf.php
Ein freies Plugin zu Acrobat muss heruntergeladen und installiert werden, um die
dortigen PDFs zu lesen.
NR: Minimum der Funktion f(x, y, . . . ) (43)Ein Algorithmus ( Downhill = Simplex = Nelder Mead ) ohne Gradient ∇f
x
y
f(x,y)
ab
c a'
Beispiel mit Mathematica
f[x_?NumberQ, y_?NumberQ] := If[x <0|| y > 0, y^2+(x-2)^2+1/(0.01+0.01*x^2 + 0.01*y^2),
y^2+(x-2)^2-1/(0.02+ 0.02 x^2 +0.01*y^2) ];
ContourPlot[f[x, y], {x, -6, 6}, {y, -4, 6}, ColorFunction -> "DarkRainbow", ContourLabels -> All]
NMinimize[{f[x, y], x < 3}, {x, y},
Method -> {"NelderMead", "InitialPoints" -> {{-6, 6} {-6, 5} {-5, 5}}}, MaxIterations -> 100]
Das Minimum hinter dem Berg bei x = y = 0wurde entdeckt!
siehe auch Bibliothek MINPACK unter Netlibhttps://en.wikipedia.org/wiki/Netlib
NR: Modeling of Data (48)Least Squares as a Maximum Likelihood Estimator
least squares fitrobust straight-line fit
y = a x + b
y = a x + b
x
y
1 2 3 4 5
10
15
20
25
x
y
(x , y )ii
y=f(x)
y=f(x)
min∑i
(yi − f(xi))2 oder min
∑i
|yi − f(xi)|
Approximation durch die Funktion y(x) = a+ bxc mit a, b, c =?
data= Table[{x, 3+2x+Random[Real,{-1, 1}]}, {x, 0.25, 5, 0.25}]; AppendTo[data, {5.25, 27.0}];
LeastSquaresError[data_, f_ ] := Sum[(f[data[[j, 1]]] - data[[j, 2]])^2, {j, 1, Length[data]}];
LeastAbsError[data_, f_ ] := Sum[Abs[f[data[[j, 1]]] - data[[j, 2]] ], {j, 1, Length[data]}];
fitting1= a+b x^c/. FindMinimum[ LeastSquaresError[data, a+b #^c &],{{a,3},{b,2},{c,1}}][[2]];
fitting2= a+b x^c/. FindMinimum[ LeastAbsError[data, a+b #^c &], {{a,3}, {b,2}, {c,1}}][[2]];
Show[ ListPlot[data, PlotStyle -> PointSize[0.025], PlotRange -> All],
Plot[fitting1, {x, -5.25, 5.25}, PlotRange -> All] ,
Plot[fitting2, {x, -5.25, 5.25}, PlotRange -> All] ]
