Lineare Algebra

ThemengebieteI Vektoren, MatrizenI Eigenwerte und EigenvektorenI Invertieren von MatrizenI Lösen linearer Gleichungssysteme

I Numerisches lösen von DGL

∆f (x) = g(x) (1)

nachdem ∆ diskretisiert wurdeI Lösen von Integralgleichungen

f (x) = g(x) +

1∫0

k(x , y)f (y)dy (2)

I Lineare Fits

VektornormSei X ein komplexer bzw. reeller Raum, und |.| eine lineareAbbildung von

|.| : X → C(R)

mitI |x | ≥ 0 (Positivität)I |x | = 0↔ x = 0 (Definitheit)I |α · x | = |α| · |x | ∀x ∈ X , ∀α ∈ C (Homogenität)I |x + y | ≤ |x |+ |y | ∀x , y ∈ X (Dreiecksungleichung)

so nennt man |.| eine Norm.Beispiel:

|x |1 =n∑

i=1

|xi | (Betragsssummennorm) (3)

|x |2 =

√√√√ n∑i=i

x2i (Euklidische Norm) (4)

SkalarproduktSeien x,y ∈ Cn(Rn) dann ist

x · y =n∑

i=1

x∗i yi im komplexen (5)

x · y =n∑

i=1

xiyi im rellen (6)

das Skalarprodukt der beiden Vektoren.In R3 gilt

x · y = |x||y| cos∠(x,y) (7)

Oder in Matrixschreibweise

〈x |y〉 =n∑

i=1

x∗1iyi1 (8)

Darstellung von Vektoren und Matrizen

Typischerweise werden Vektoren und Matrizen als Listen bzw.verschachtelte Listen dargestellt

I Python:a=[[1,2],[2,0]]

I C (z.B. GSL):m->data[0*m->tda+0]=1.;gsl_matrix_set(m,0,1,2.);gsl_matrix_set(m,1,0,2.);gsl_matrix_set(m,1,1,0.);

I Mathematica:a={ { 1,2},{2,0} }

Inversion der Matrix in Python (numpy)/MMAPython

import numpy as np

a=[[ 1,2],[2,0]]ainv=np.linalg.inv(a)

print(ainv)print(np.dot(ainv,a))

Mathematica

a={ {1,2}, {2,0} }ainv=Inverse[a]ainv.a

Output

MathematicaScript -print all -f ./testmath{{ 1, 2}, {2, 0}}{{ 0, 1/2}, {1/2, -1/4}}{{ 1, 0}, {0, 1}}

Inversion der Matrix in gsl#include <stdio.h>#include ‘‘gsl/gsl_matrix.h’’#include ‘‘gsl/gsl_linalg.h’’

int main(){

gsl_matrix *m;gsl_matrix *minv;gsl_permutation *p;int s;

p=gsl_permutation_calloc(2);m=gsl_matrix_alloc(2,2);minv=gsl_matrix_alloc(2,2);

m->data[0*m->tda+0]=1.;gsl_matrix_set(m,0,1,2.);gsl_matrix_set(m,1,0,2.);gsl_matrix_set(m,1,1,0.);

gsl_linalg_LU_decomp (m, p, &s);gsl_linalg_LU_invert(m,p,minv);

printf(‘‘The inverse is\n’’);

printf(‘‘%lf , %lf\n’’,gsl_matrix_get(minv,0,0),gsl_matrix_get(minv,0,1));printf(‘‘%lf , %lf\n’’,gsl_matrix_get(minv,1,0),gsl_matrix_get(minv,1,1));

return 0;

}

Vektoren und Matrizen in der Physik

Praktisch werden in der Physik ständig Vektoren und Matrizengebraucht:

I Mechanik (Kräftevektoren)I Geometrische OptikI Quantenmechanik

Einige Definitionen sind hierbei gebräuchlich, insbesonderehinsichtlich Darstellung von Matrizen über ihre Einträge.

MatrizenDarstellung einer Matrix über Einträge (Koordinatenform)

An×m = aij mit i ≤ n ,j ≤ m (9)

In dieser Schreibweise ist das Produkt zweier Matrizen:

A · B =∑

j

aijbjm = aijbjm = cim (10)

Damit sind auch die beim Skalarprodukt eingegführten brasund kets

〈x |i = x1i (Zeilenvektor) (11)|y〉i = yi1 (Spaltenvektor) (12)

Man lässt häufig die Summenzeichen weg und meint das übergleiche Indices summiert wird ≡ EinsteinscheSummenkonvention.

Einige Operationen mit Matrizen:

I Transposition

AT = (aij)T = aji

(A · B)T = (aijbjm)T = bmjaji = BT · AT (13)

I Adjungierte

(AT )∗ = A† =(

(aTij )∗)

= a∗ji (14)

I Spur

Tr(A) = aii (15)

I Determinante

detA = εi1...ina1i1 . . . anin (16)

I ε ist der Levi-Civita-Tensor mit den Eigenschaften

εijk ... =

+1 falls (i , j , k , . . . ) gerade Permutation−1 falls (i , j , k , . . . ) ungerade Permutation0 wenn mindesten zwei Indices gleich sind

(17)

I Häufig auch Kreuzprodukt zweier Vektoren inMatrixschreibweise

(~a× ~b)i = εijkajbk (18)

Hier das Bsp.

~a · (~a× ~b) = aiεijkajbk = −εjikaiajbk = 0 (19)

Besondere Matrizen, z.B.I Hermitesch (Symmetrisch)

A† = A (20)

(aij)† = a∗ji = aij (21)

I Unitär (Orthogonal)

A† = A−1 (22)

Weshalb sind diese beiden Eigenschaften so besonders?I Hermitesch⇒ reelle Eigenwerte (physik. Observable)I Unitär⇒ Normerhaltende Transformationen

(Basiswechsel), alle EW liegen auf Einheitskreis

Zeige beide Eigenschaften in HilbertraumHermitesch (sei A Eigenzustand zu O mit Eigenwert a), d.h.

O|A〉 = a|A〉

Zunächst versteht sich in bra-ket-Schreibweise

O|A〉† =(

OijAj1

)†=(

A∗1jO∗ji

)= 〈A|O† (23)

〈A|O|A〉 = a〈A|A〉〈A|O|A〉 = 〈A|O†|A〉∗ = 〈A|O|A〉∗

= a∗〈A|A〉⇒ a = a∗ (24)

Zeige beide Eigenschaften in Hilbertraum

Operator U ist unitär.

〈A|A〉 = 〈A|U†U︸︷︷︸1

|A〉

=(〈A|U†

)︸ ︷︷ ︸

u∗〈A|

(U|A〉

)︸ ︷︷ ︸

u|A〉

= u∗u〈A|A〉⇒ u∗u = 1 (25)

Eigenwerte haben Norm 1⇒ u = eiφ

Häufige Frage: Finde eine unitäre Transformation, sodass eineMatrix Diagonalform annimmt (z.B. Hauptachsentrafo)

Eigenwerte und EigenvektorenI Eigenwertgleichung für Operator O mit EV Ai und EW ni

OAi = niAi (26)

I Bestimmung über charakteristisches Polynom

det(O − λI) = 0 (27)

Nullstellen sind die Eigenwerte.I Trafo

UOU† = O (28)

wobei U zeilenweise EV sind

U =

A1,1 A1,2 . . . A1,n...

......

...An,1 An,1 . . . An,n

(29)

I Maximaler EW (betragsmässig) wird spektraler Radiusgenannt

ρ(A) = max{|λ||λ ∈ σ(A)} (30)

Beispiel:

M =

(2 11 −2

)det(M − λI) = λ2 − 5 = 0

λ1,2 = ±√

5(M − λ1,2I)v1,2 = 0

v1 = (2−√

5,1)

v2 = (2 +√

5,1)

U =

(2−√

5 12 +√

5 1

)UMU† =

(−10

(−2 +

√5)

00 10

(2 +√

5) ) (31)

MatrixzerlegungenI Für alle Matrizen A ∈ Cn×n existiert eine unitäre Matrix

U ∈ Cn×n derart, dass

M = UAU† (32)

eine rechte obere Dreiecksmatrix ist.I Sollte A hermitesch sein ist M diagonal.I Definiere die Koniditionszahl einer Matrix A

cond(A) = |A|a|A−1|a (33)

mit |A|a := sup|x |a=1

|Ax |a induzierte Matrixnorm. Je größer

diese Zahl, desto schlechter konvergieren iterative LöserI Schreibe Matrix A ∈ Rn×n um

UAV T = diag{σ1, . . . , σn} (34)

mit U und V orthogonale Matrizen und

0 < σ1 ≤ · · · ≤ σn (Singulärwerte) (35)

SingulärwertzerlegungI Bestimmungsgleichung der Singulärwerte

det(AT A− σ2i I) = 0 (36)

I SW sind EW von AT AI Damit ist V unitäre (orthogonale) Trafo

VAT AV T = diag{σ21, . . . , σ

2n} (37)

zeilenweise EV von AT AI U lässt sich über das Inverse bestimmen

U = diag{σ1, . . . , σn}(

AV T)−1

(38)

I Konditionszahl ist dann

cond(A) =σn

σ1(39)

Singulärwertzerlegung hat viele AnwendungenI Einfache Bildung der inversen, da U und V unitär

(orthogonal sind) braucht nur noch die Kehrwerte derDiagonalmatrix

I Bestimmung der PseudoinversenI Lösung des LeastSquare-Problems

Es existieren noch weitere Zerlegungen (LR, QR, Choleskyusw.)

Inversion von Matrizen

Generell gehen wir von quadratischen Matrizen aus.I Direkte Verfahren z.B. über Adjunkte, SVD oder andere

ZerlegungenI Iterative Verfahren (Splitting-Verfahren z.B. Jacobi)

Frage der Eindeutigkeit von Lösungen zu Ax = bI Inverse existiert wenn det(A) 6= 0

I...

A−1 =adj(A)

det(A)(40)

Verfahren über Adjunkte

Die Adjunkte ist definiert

adj(A) =

a11 a12 . . . a1na21 a22 . . . a2n

.... . .

an1 an2 . . . ann

T

aij = (−1)i+jdet

a11 . . . a1j−1 a1j+1 . . . a1n...

. . ....

......

ai−11 . . . ai−1j−1 ai−1j+1 . . . ai−1nai+11 . . . ai+1j−1 ai+1j+1 . . . ai+1n

.... . .

......

. . .an1 . . . anj−1 anj+1 . . . ann

(41)

Python

import numpy as np

def adjunkte(a):b=np.zeros( (len(a),len(a)))for i in range(len(a)):

for j in range(len(a)):subdet=np.linalg.det(np.delete(np.delete(a,i,0),j,1))b[i,j]=(-1)**(i+j)*subdet

return np.transpose(b)

def inv(a):return adjunkte(a)/np.linalg.det(a)

Über SVD

I Nutze aus, dass in der Zerlegung

USV T = A (42)

U und T unitär (orhtogonal) sind

A−1 =(

V T)−1

S−1U−1

= VS−1UT (43)

I Die inverse der Diagonalmatrix sind einfach die Kehrwertein der Diagonalen

Splitting-VerfahrenZiel ist es eine Fixpunktiteration zu bekommen

xn+1 = F (xn) (44)

Wir lösen allgemein:

Ax = b (45)

I Iterationsverfahren ist geg. durch AbbildungF : Cn × Cn → Cn und heißt linear wenn

F (x,b) = Mx + Nb (46)

I x ∈ Cn bezeichnen wir als Fixpunkt der Abb.F : Cn × Cn → Cn, wenn

x = F (x,b) (47)

mit b ∈ Cn.I F heißt konsistent zur Matrix A, wenn

x = limn→∞

xn = limn→∞

F (xn−1,b) unabh. vom Startwert

I Iterationsverfahren genau dann konvergent wenn

ρ(M) < 1 (48)

der Spektralradius der Iterationsmatrix < 1 ist.I Splitting-Methode

A = B + (A− B)

⇒ Ax = bBx = (B − A)x + b

→ x = B−1(B − A)︸ ︷︷ ︸M

x + B−1︸︷︷︸N

b (49)

I Lineares Iterationsverfahren

xn+1 = F (xn,b) = Mxn + Nb (50)

Jacobi-Verfahren

I Wenn A nichtverschwindende Diagonalelemente aufweistaii 6= 0 benutze

D = diag{a11, . . . ,ann}M = D−1(D − A)

N = D−1

Iteration

x = D−1(D − A)x + D−1b (51)

I Wiederum Inversion der Diagonalmatrix ist sehr einfach.I Konvergenz wenn spektraler Radius von M < 1

Weitere Verfahren

I Gauß-Seidel ≡ SplittingverfahrenI MehrgitteverfahrenI Projektive Verfahren (Krylov-Unterraum-Methoden)

Nichtquadratische Matrizen

Was macht man bei Nichtquadratischen Matrizen A ∈ Rm×n?I Erzeuge quadratische Matrix

Ax = b

⇒ AT Ax = AT b

x =(

AT A)−1

AT︸ ︷︷ ︸P

b (52)

P ist pseudoinverseI Sollte AT A nicht invertierbar sein definiert man die

Pseudoinverse über die SVD bei der man die SW < εvernachlässigt.

Least-Squares-Problem

I Beziehung zum Least-Squares-Problem

r = Ax − y (53)

definiert das Residuum, dann ist das Abstandsquadrat

|r |2 = rT r = (xT AT − yT )(Ax − y)

= xT AT Ax − xT AT y − yT Ax + yT y

= xT AT Ax − 2yT Ax + yT y (54)

I Finde das Minimum (1.Abl Null setzen)

∇|r |2 = 2AT Ax − 2yAT = 0

⇒ AT Ax = AT y (55)