Numerische Behandlung des Eigenwertproblems

Numerische Behandlung des Eigenwertproblems

Zusammenfassung

Das Ziel dieses Vortrages ist, zwei gute Methoden für die numerische Bestimmung der Eigenwerte zu zeigen und wie man diese mit Matlab anwenden kann (implementieren).

Namentlich wären das die Vektoriteration und die QR-Zerlegung.

0. Einleitung

Wie wir bereits gesehen haben, benötigt man zum Lösen von den verschiedensten Problemen oft die Eigenwerte (EW) einer Matrix; sei es um die Hauptspannungen einer belasteten Probe auszurechnen oder eine Matrix in der Normalform zu schreiben.

Deshalb möchte man die Eigenwerte möglichst gut und einfach berechnen können. In der Praxis geschieht dies natürlich numerisch mittels Computern.

Mit der Methode des charakteristischen Polynoms, die wir ja meist anwenden, hat der Computer seine Mühen, besonders wenn die Eigenwerte sehr nahe beieinander liegen. In diesem Fall können unverhältnismässig grosse Fehler entstehen, was gelinde gesagt, nicht wünschenswert ist.

Ein Auszug vom Buch „Lineare Algebra“ von K. Nipp auf der Seite 215 soll diese Problematik verdeutlichen:

9598 . 3 6083 . 7 8142 . 4

) 2 )(

2 )(

4 . 1 ( ) (

2 3

3

+

− +

−

=

−

−

−

−

=

λ λ

λ

λ λ

λ λ

P

Das Polynom mit 5-stelliger Numerik ausgerechnet, ist mit 3-stelligen Fehlern behaftet.

15

* 3 3 3

3

* 2 2 2

3

* 1 1 1

10

* 3 . 1 ] [

10

* 6 . 1 ] [

10

* 6 . 1 ] [

−

−

−

≈

−

=

∆

≈

−

=

∆

≈

−

=

∆

λ λ

λ λ

λ λ

Wobei das die exakten Eigenwerte sind

2 2

4 .

1 _{2 3}

1 = λ = λ =

λ

Zusammenfassend kann man sagen, dass die Methode des charakteristischen Polynoms folgende Nachteile besitzt:

- Bei Polynomen höheren Grades ist die Bestimmung der Nullstellen schwer und fehleranfällig

- imaginäre Lösungen können auftreten

- relativ grosse Fehler bei nah zusammenliegenden EWs - (es werden alle EW ausgerechnet)

Es sind also andere Methoden von Nöten, um die EW zu bestimmen.

1. Vektoriteration

1.1.1 Vektoriteration für symmetrische Matrizen

Die Vektoriteration ist ein Verfahren zur Berechung von Eigenvektoren und Eigenwerten. Es wird dabei jeweils nur der betragsmässig kleinste Eigenwert berechnet, in der Praxis reicht das aber oft aus.

Wir betrachten hier nur die Vektoriteration für symmetrische Matrizen, am Beispiel einer 2x2 Matrix A. Für eine solche Matrix A gilt:

− Alle Eigenwerte sind reell

− Es existiert eine orthonormale Eigenbasis zu A

Aufgrund dieser Voraussetzungen lässt sich ein beliebiger 2x1 Vektor x^a als

Linearkombination der beiden normierten Eigenvektoren x⁽¹⁾ und x⁽²⁾ der Matrix A schreiben:

x^a = c₁^a x⁽¹⁾ + c₂^a x⁽²⁾

Es gelten ausserdem folgende Identitäten:

A x⁽¹⁾ = λ₁ x⁽¹⁾

A x⁽²⁾ = λ₂ x⁽²⁾ (Definition des Eigenwerts) A x^a = A c1a

x⁽¹⁾ + A c2a

x⁽²⁾ = λ₁ c1a

x⁽¹⁾ + λ₂ c2a

x⁽²⁾ 1.1.2 erste Iteration

Wir definieren nun eine Folge von Vektoren x⁰, x¹, x², … , x^k.

Der obere Index ohne Klammer bezeichnet den jeweiligen Iterationsschritt. x⁰ ist ein zufällig gewählter Startvektor.

Es soll nun gelten:

A x¹ = x⁰ A x² = x¹ A x³ = x² Usw.

Dies ist unsere Iterationsvorschrift. Für die Iteration muss in jedem Schritt eine Gleichung A x^k = x^k-1

mit Hilfe einer LR-Zerlegung gelöst werden, dies führt am Ende zu einem Vektor x^k. 1.1.3 Eigenschaften dieser Iteration

Unter Anwendung obiger Identitäten lässt sich die Gleichung schreiben als:

λ₁ c₁¹ x⁽¹⁾ + λ₂ c₂¹ x⁽²⁾ = c₁⁰ x⁽¹⁾ + c₂⁰ x⁽²⁾ (λ1 c11

- c10

) x⁽¹⁾ + (λ2 c21

– c20

) x⁽²⁾ = 0

Da x⁽¹⁾ und x⁽²⁾ die Eigenbasis von A bilden, sind sie linear unabhängig und die obige Gleichung wird nur Null, wenn die Koeffizienten gleich Null sind:

(λ₁ c₁¹ - c₁⁰) = 0 (λ₂ c₂¹ - c₂⁰) = 0

⇒ c11

= (1/ λ1) c10

c21

= (1/ λ2) c20

Der Iterationsvorschrift zufolge lässt sich der Vektor:

x¹ = c11

x⁽¹⁾ + c21

x⁽²⁾ also schreiben als:

x¹ = (1/ λ₁) c₁⁰ x⁽¹⁾ + (1/ λ₂) c₂⁰ x⁽²⁾

In einem 2. Iterationsschritt können wir x² bestimmen:

A x² = x¹

⇒ λ1 c12

x⁽¹⁾ + λ2 c22

x⁽²⁾ = c11

x⁽¹⁾ + c21

x⁽²⁾

Mit der gleichen Rechung wie oben lassen sich wieder c₁² und c₂² durch c₁¹ und c₂¹ ausdrücken:

c12

= (1/ λ₁) c11

c₂² = (1/ λ₂) c₂¹

⇒ c₁² = (1/ λ₁)² c₁⁰ c22

= (1/ λ2)² c20

Daraus folgt für den Vektor x²: x² = (1/ λ1)² c10

x⁽¹⁾ + (1/ λ2)² c20

x⁽²⁾ = (1/ λ1)² (c10

x⁽¹⁾ + (λ1/ λ2)² c20

x⁽²⁾) Für k Iterationsschritte gilt:

x^k = (1/ λ1)^k (c10

x⁽¹⁾ + (λ1/ λ2)^k c20

x⁽²⁾)

Wir definieren nun λ₁ als den betragsmässig kleinsten Eigenvektor der Matrix A.

Dann konvergiert der zweite Term der Summe für k → ∞ gegen Null, da (λ1/ λ2) < 1.

Es gilt dann:

x^k = (1/ λ₁)^k c₁⁰ x⁽¹⁾ = d x⁽¹⁾

d ist dabei nur ein Skalar, d.h. die Iteration führt zu einem x^k , welches gegen einen Vektor in Richtung des Eigenvektors konvergiert.

1.1.4 zweite Iteration

Mit diesem Wissen modifizieren wir die Iterationsvorschrift so, dass wir jeden iterierten Vektor normieren. Dies spielt vor allem für die spätere Berechnung von λ1 eine Rolle. Wir starten mit einem normierten x⁰:

A x¹ = x⁰

A x² = (x¹ / ||x¹||) A x³ = (x² / ||x²||) usw.

Für k → ∞ gilt dann:

(x^k/ ||x^k||)= x⁽¹⁾

Somit kann man mit dieser Iteration den Eigenvektor zum kleinsten Eigenwert der Matrix A berechnen. Die Iteration wird solange durchgeführt bis die Vektoren (x^k/ ||x^k||) und (x^k-1/||x^k-

1||) bis zur gewünschten Genauigkeit übereinstimmen.

Für die Berechnung des Eigenwertes λ₁ zu diesem Eigenvektor, betrachten wir die folgende Gleichung:

Für k gross genug gilt:

(x^k/ ||x^k||)= x⁽¹⁾ A(x^k/ ||x^k||)= Ax⁽¹⁾ (Ax^k)/ ||x^k||= λ1x⁽¹⁾ (x^k-1/ ||x^k-1||)/||x^k|| = λ₁x⁽¹⁾

||((x^k-1/ ||x^k-1||)/||x^k||)|| = ||λ₁x⁽¹⁾||

1/||x^k|| = λ1 (Da x⁽¹⁾ und (x^k-1/ ||x^k-1||) die Norm 1 besitzen)

Somit kann man durch die Vektoriteration bis zum k-ten Schritt den kleinsten Eigenwert und den dazugehörigen Eigenvektor berechnen.

1.2 Die Implementierung der Vektoriteration in Matlab 1.2.1 Vorgehensweise

Um ein Algorithmus erfolgreich zu implementieren, ist es sicher sinnvoll die einzelnen Schritte des Verfahrens in der logischen Reihenfolge zu skizzieren:

Zur Erinnerung:

Gegeben ist: A sei eine beliebige, symmetrische ( A = A^T ) Matrix vom Rang r und x0 sei ein beliebiger Spaltenvektor ( r x 1), jedoch kein EV von A Gesucht ist: Betragsmässig kleinster EW von A

Input: Dimension Toleranz Matrix

x₀ normieren

Gleichungssystem Ax = x₀ nach x lösen

x normieren

x mit x₀ vergleichen Genügend genau?

Nein

Ja Output: kl. Eigenwert

Eigenvektor

1.2.2 Das Programm

Das Schema aus 1.2.1 in die Sprache von Matlab übersetzt wäre dann:

% Vektoriteration Titel

d = 3 legt die Dimension fest

A = randn(d,d); erstellt eine dxd Matrix, hier 3x3

A = A'*A; macht sie symmetrisch

A = round(10*A) rundet Werte (zwischen 0 und 10)

x0 = randn(d,1); erstellt zufälliger Spaltenvektor

x0 = round(x0) rundet den Spaltenvektor

x0 = x0/norm(x0); x0 wird durch seine Norm ersetzt (Normierung)

for i=1:100 Begin der Schleife (*)

x = A\x0; neuer Vektor x

n = norm(x);

x = x/n; normiert den Vektor x

if (abs(x-x0)<1.e-5) Abbruchsbedingung, Toleranz break;

else

x0 = x; Definiert x als neues x0

end (es geht zurück zu (*))

end

fprintf('Anzahl von Schritten %i\n',i) Ausgabe der benötigten Schritte fprintf('Eigenwert ist %f\n',1/n) zeigt den EW

x zeigt den EV

[V D] = eig(A); Kontrollrechnung

D(1,1) zeigt den EW

V(:,1) zeigt den EV

Es muss nicht jeder einzelne Tag (Code) verstanden werden, sondern die Vorgehensweise ist entscheidend.

Zeichenerklärung:

% kommentiert eine Zeile aus (deaktiviert)

; Unterdrückt die Anzeige dieser Zeile im Output fprintf(' ') Textfeld

for Erstellt i Schleifen (i geht von 1 bis 100)

break Stoppt die for-Schleife

if else Wenn (if) die Bedingung erfüllt ist, dann mache dies, sonst (else) mach das.

randn Erzeugt eine Zufallszahl

norm() Normiert einen Vektor

round() Rundet Zahlen

[V D] = eig(A) Erstellt eine Matrix V bestehend aus den EV und einer Diagonalmatrix D mit den EW.

Hier wurde eine Zufallsmatrix der Dimension d erstellt. Man kann auch eine selbstgewählte (symmetrische) Matrix verwenden.

Beispielsweise: A =[ 1 2 ; 2 1]

Die Toleranz oder die Mindestgenauigkeit wird bei if mit: (abs(x-x0)<1.e-5) bestimmt. Das Programm läuft bis diese Genauigkeit erreicht wurde (oder i=100).

Dass das Programm tatsächlich funktioniert, zeigt das Beispiel im Anhang 3.1

1.2.3 Vergleich: Methode des charakteristischen Polynoms und der Vektoriteration

Wir kommen zurück zur Methode des charakteristischen Polynoms und möchten diese mit der Vektoriteration vergleichen. Dazu wurde auch die Methode des charakteristischen Polynoms in Matlab implementiert. (Der Quellcode befindet sich im Anhang 3.2)

Da es schwer ist Rechenschritt an sich zu zählen, möchten wir die benötigte Zeit der beiden Algorithmen vergleichen. Die Zeit ist schlussendlich ja das entscheidende Effizienzkriterium.

Dazu wurde die Zeit mit der Funktion ’tic’ ’toc’ vom Computer selbst gemessen. Natürlich mussten die beiden Algorithmen die gleiche Zufallsmatrix bearbeiten.

Es wurden verschiedene d für die Dimension verwendet. Dabei hat sich herauskristallisiert, dass bei kleinen d beide Algorithmen praktisch gleich schnell sind, hingegen bei grösseren d, die Vektoriteration deutlich schneller ist.

Jedoch muss man eingestehen, dass bei der Methode des charakteristischen Polynoms alle Eigenwerte ausgerechnet wurden.

1.2.4 Gesamtbeurteilung der Vektoriteration

Somit hat die Vektoriteration gegenüber der Methode des charakteristischen Polynoms folgende Vorteile:

- Liefert beliebig genaue Resultate

- Geringe Fehleranfälligkeit, keine komplexen Zahlen - Hohe Effizienz

- Es wird nur das berechnet, was man wirklich braucht.

Eine Schwachstelle des Algorithmus ist, wenn der zufällige Vektor gerade ein Eigenvektor der Matrix ist. Dies ist aber unwahrscheinlich und deshalb nicht weiter störend.

Da es aber auch vielfach Fälle gibt, wo man alle Eigenvektoren und alle Eigenwerte einer Matrix ausrechnen will, braucht man noch einen zweiten Algorithmus.

2. QR Algorithmus 2.1 Ziel QR-Algorithmus:

Finden alle Eigenwerte und Eigenvektoren ( auch Komplexe) einer allgemeinen reellen n x n- Matrix.

2.1.2 Vorwissen aus QR-Zerlegung:

Zu jeder m x n-Matrix A, mit m > n, existiert eine orthogonale m x m-Matrix Q, so dass gilt:

A = QR mit R = ( R₀/ 0)

Wobei R0 eine n x n-Rechtsdreiecksmatrix ist und 0 die (m-n) x n-Nullmatrix.

2.1.3 Theorie:

Zwei quadratischen Matrizen A und B werden ähnlich genannt, falls es eine reguläre n x n- Matrix T gibt, so dass:

B = T^-1AT

In diesen Fall, haben A und B die gleiche Eigenwerte, nämlich, wenn λ Eigenwerte von A ist, und wenn x seine Eigenvektor ist, es gilt:

BT^-1x = T^-1Ax = λ T^-1x Beweis:

Ist λ Eigenwert von A und x seine Eigenvektor, es gilt Ax1 = λx1

Wir müssen beweisen, dass es gilt:

B x2 = λx2 x1 ≠ x2

Wir wissen auch dass es gilt: B = T^-1AT

Wir können probieren, diese Gleichung zu benutzen, aber es folgt:

T-1AT x = λ x !!!!! Diese ist falsch!!! Das T zwischen A und x muss verschwinden!

Um diese T zu eliminieren, können wir B multiplizieren nach T^-1. B T^-1 = T^-1AT T^-1 = T^-1AI = T^-1A

Jetzt können wir T^-1A benutzen:

T^-1Ax = λ T^-1x = BT^-1x

Wir sehen, dass die Eigenwert von A und B ist gleich, aber die Eigenvektoren sind verschiedene:

Eigenvektor von A = x= x1

Eigenvektor von B = T^-1x = x2

Also λ ist auch Eigenwerte von B und seine Eigenvektor ist T^-1x.

Die Methoden, die sie gleichzeitig alle Eigenwerte von einer Matrix zu aproximieren erlauben, sind generell auf die Idee begründet sich die Anfangsmatrix in eine ähnliche

Dreiecksmatrix zu verwandeln. Davon entsteht, dass die Eigenwerte von die Diagonalelement dargestellt werden.

Die QR-Algorithmus ist das mehr benutzte Methode, um die Eigenwerte zu berechnen, auch, weil das wirksamste ist.

2.1.4 QR-Algorithmus:

Zerlege: A1 = Q1R1

Bilde A₂ = R₁Q₁ [ = Q₁^TAQ₁]

Dieser zwei Schritte sind die Basis des QR-Algorithmus, findet einmal A1, es muss dem gleichen Verfahren folgen, um A₃ A₄,… A_n zu finden.

Je höher n wird, desto genauer sind die gefundene Eigenwerte.

In allgemein für das QR-Algorithmus gilt:

Ak = QkRk

und

Ak+1 = RkQk

Von diese zwei Beziehungen folgt, dass:

Rk= Qk-1

Ak

A_k+1 = Q_k^-1A_kQ_k =Q_k^TA_kQ_k

Also die Matrizen der Folge {Ak} sind alle ähnliche, und also haben die gleiche Eigenvektor.

Die Matrix Ak für k ∞ konvergiert zu eine Dreiecksmatrix.

2.1.5 Die Eigenvektoren:

Um die Eigenvektor x^(j) der Matrix A zu finden, muss man die folgende Gleichung zu lösen:

x^(j) = Q y^(j)

Wobei y^(j) ist die Eigenvektor der Matrix A_k (und es ist einfach zu berechnen!) und Q = Q1*Q2*…..*Qk ist die Produkt aller orthogonale Matrizen Qi, die wir benutzen haben, um das QR-Algorithmus zu führen.

2.1.6 Einfach Beispiel: (mit MATLAB durchgeführt)

Bestimmen alle Eigenwerte dieser Matrix mit das QR-Algorithmus:

A₁ = [ 3 4 ; 4 1]

Zerlegung von A1:

A1 = Q1R1 = [ -0.6000 -0.8000 ; -0.8000 -0.6000] * [ -5.0000 -3.2000 ; 0 -2.6000]

Bildung von A2:

A2 = R1Q1 = [5.5600 2.0800 ; 2.0800 -1.5600]

Zerlegung von A2:

A2 = Q2R2 =[-0.9366 -0.3504 ; -0.3504 0.9366] * [-5.9363 -1.4015 ; 0 -2.1899]

Bildung von A₃:

A3=R2Q2= [6.0511 0.7673 ; 0.7673 -2.0511]

Zerlegung von A₃:

A3=Q3R3= [-0.9921 -0.1258 ; -0.1258 0.9921] *[-6.0996 -0.5032 ; 0 -2.1313]

Bildung von A₄:

A₄= R₃Q₃= [6.1144 0.2681 ; 0.2681 -2.1144]

….

….

Bildung von A7:

A₇= R₆Q₆ = [6.1231 0.0112 ; 0.0112 -2.1231]

Jetzt können wir stoppen, weil wir eine hinlänglich Genauigkeit erreichen haben.

Die gefundene Eigenwerte sind deshalb 6.1231 und -2.1231.

Wenn wir mit der Methode der Determinante die Eigenwerte berechnen, erhalten wir die gleiche Resultat.

Berechnen nun die Eigenvektoren:

Q = [0.7891 -0.6143 ; 0.9295 -0.3688]

y⁽¹⁾ = [ 0.0014 ; -1]

y⁽²⁾= [ -1 ; -0.0014]

Benutzen wir jetzt die Formel x^(j) = Q y^(j) : x⁽¹⁾ = [ 0.6154 ; 0.3701]

x⁽²⁾ = [ -0.7882 ; -0.9290]

2.2 Beispiele zum QR-Algorithmus [mit Matlab]

2.2.1 Einfache 2x2-Matrix

A = [1 2;2 1]

A = 1 2 2 1

Kontrolle eig(A)

ans = -1 3

Eigenwerte: -1 und 3

2.2.1.1 Schrittweise Berechnung

[Q R] = qr(A) Q =

-0.4472 -0.8944 -0.8944 0.4472

R = -2.2361 -1.7889 0 -1.3416 A= Q*R

A = 1.0000 2.0000 2.0000 1.0000

A1= R*Q

A1 = 2.6000 1.2000 1.2000 -0.6000

[Q1 R1]= qr(A1) Q1 = -0.9080 -0.4191

-0.4191 0.9080 R1 = -2.8636 -0.8381

0 -1.0476

A3=R2*Q2

A3 = 2.9945 0.1479 0.1479 -0.9945

[Q3 R3]= qr(A3) Q3 = -0.9988 -0.0493 -0.0493 0.9988 R3 = -2.9982 -0.0987

0 -1.0006

A2= R1*Q1

A2 = 2.9512 0.4390 0.4390 -0.9512

[Q2 R2]=qr(A2) Q2 = -0.9891 -0.1471

-0.1471 0.9891 R2 = -2.9837 -0.2943

0 -1.0055

A4=R3*Q3

A4 = 2.9994 0.0494 0.0494 -0.9994 usw.

ev = 3.0000 -1.0000 i = 12

2.2.1.2 Direkte Berechnung mit Matlab

A = [1 2;2 1]

Matrix A = 1 2 2 1

Genauigkeit oder Toleranz tol=1.e-5;

Maximaler Lauf maxn=100;

Titel

fprintf('QR Algorithmus\n') QR Algorithmus

for i=1:maxn [Q R] = qr(A);

A= Q'*A*Q;

Abbruchbedingung if (abs(A(2,1))<tol) Eigenwerte von A ev= diag(A);

break

end

end

2.2.2 Beliebige 5x5 Matrix

tol= 1.e-5;

maxn=700;

Beliebige 5x5-Matrix B= randn(5,5);

Symmetrische Matrix B=B'*B

Matrix B = 4.8707 0.3102 -2.5578 0.3060 -2.2967 0.3102 1.5799 -1.2100 -0.9150 -0.0056 -2.5578 -1.2100 4.9952 0.6148 -0.4092 0.3060 -0.9150 0.6148 2.1305 -0.9062 -2.2967 -0.0056 -0.4092 -0.9062 2.7951

Titel

fprintf('QR Algorithmus\n') QR Algorithmus

for i=1:maxn [Q R]=qr(B);

B=Q'*B*Q;

if (abs(B(2,1))<tol)

Kontrolle eig(B)

ans = 8.0322 5.0220 1.9916 1.0449 0.2808 Eigenwerte von B

ev=diag(B);

break

end

end

i = 26

ev = 8.0322 5.0220 1.9916 1.0449 0.2808

3. Anhang

3.1 Matlab Ausgabe einer Vektoriteration Ein Beispiel einer Vektoriteration mit Matlab

d = Gewählte Dimension

4

A = Zufällige symmetrische Matrix

80 -5 -52 -51 -5 9 -4 13 -52 -4 66 36 -51 13 36 50

x0 = Zufälliger Startvektor

0 1 -1 -1

x = Erste Iteration

-0.1389 0.8159 0.2159 -0.5181

x = Zweite Iteration

-0.1264 0.8170 0.2312 -0.5129

x = Dritte Iteration

-0.1259 0.8171 0.2316 -0.5126

x = Vierte Iteration

-0.1259 0.8172 0.2316 -0.5126

x = Fünfte Iteration

-0.1259 0.8172 0.2316 -0.5126

Anzahl von Schritten 5

Eigenwert ist 0.481955 Eigenwert (Lösung)

x = Eigenvektor (Lösung)

-0.1259 0.8172 0.2316 -0.5126

3.2 Programm: Algorithmen-Vergleich

Dieses Programm vergleicht die Laufdauer der beiden Methoden.

% Parameter festlegen

tol = 1.e-6; % Toleranz

maxn = 100; % Maximale Anzahl Schleifen

d = 4; % Dimension

% Matrix wird erstellt

A = randn(d,d); % Erstellt eine beliebige Matrix

A = A'*A; % Wird symmetrisch gemacht

A = round(10*A)

% Eigenwerte ausrechnen % Kontrollrechnung

fprintf('Genaue Eigenwerte\n') [EV EW] = eig(A);

ev = diag(EW) % Gibt die Ews aus

% Errechnet die EW mit dem charakteristischem Polynom fprintf('Löst mit der Methode des charakteristischen Polynoms\n')

tic % Startet Zeitmessung

syms lambda

p = det(A-lambda*eye(d)); % Das charakteristische Polynom

ev = solve(p); % Findet Nullstellen

ev = eval(ev) % Gibt die Ews aus

toc % Stoppt Zeitmessung

% Vektoriteration

fprintf('Löst mit der Methode der Vektoriteration\n')

tic % Startet Zeitmessung

y0 = randn(d,1); % Zufälliger Startvektor

y0 = y0/norm(y0); % Normiert

for i=1:maxn

y = A\y0; % Nächster Vektor

n = norm(y); % Normiert

y = y/n;

if (abs(y-y0)< tol) % Prüft das Stoppkriterium

ev = 1/n;

break; % Stoppt

else

y0 = y; % Neue Itineration

end end

ev % Gibt den EW aus

toc % Stoppt Zeitmessung

return