Householder-Verfahren < Eigenwertprobleme < Numerik < Hochschule < Mathe < Vorhilfe
|
Ich habe diese Frage in keinem Forum auf anderen Internetseiten gestellt.
Ich suche Quellcode des Householder-Verfahrens (Eigenwerte + Eigen-
vektoren einer) - habe bereits 4 Varianten getestet, aber im Falle
einer 10x10-Matrix versagen die Algorithmen.
Konkret:
Der Quellcode aus dem Buch "Numerical Recipes in C - 2. Auflage"
liefert noch die besten Ergebnisse (fast alle Eigenwerte korrekt,
aber die Eigenvektoren sind falsch - habe es mit Mathematica über-
prüft). Da Mathematica keine Beanstandungen hatte, ist somit auch
von einer korrekten (geeigneten) Matrix auszugehen.
Tippfehler im Quellcode kann ich inzwischen auch ausschliessen.
THEO
|
|
| |
|
Schau mal diese Folie an. Villeicht hilft dir weiter, selber habe ich es nicht programmiert.
Quellcode habe ich aus meine Vorlesungsscript
function [A] = QR_HOUSE(A);
% input: A - (n x m) matrix
% output: A - matrix containing the necessary information of the Householder
% vectors v in the lower triangle and R in the upper triangle
[n,m] = size(A);
for k = 1:min(n-1,m)
v(k:n,1) = HOUSEHOLDER(A(k:n,k));
function [A] = HOUSEHOLDER_MULT(A,v);
% input: A - matrix
% v - Householder vector
% output: A - transformed matrix
% A = A - 2
v’v(vv’)A
beta = -2/(v’*v);
w = v’*A; % w is a line vector
A = A + beta*v*w;
function [Q] = Q_HOUSE(A);
% input: A - (n x m) matrix, QR - decomposition containing the
% necessary information of the Householder vectors v
% in the lower triangle and R in the upper triangle
% output Q - orthonormal matrix Q = [mm] Q_n-1* [/mm] .... [mm] *Q_1
[/mm]
% with [mm] Q_k [/mm] = I-2/(v’*v)*(v*v’)
n = size(A,1);
Q = eye(n);
for k = 1:min(n-1,m)
v = ones(n+1-k,1);
v(2:n+1-k) = A(k+1:n,k);
Qk = eye(n);
Qk(k:n,k:n) = eye(n+1-k) - (2/(v’*v))*(v*v’);
Q = Qk*Q;
end
A(k:n,k:m) = HOUSEHOLDER_MULT(A(k:n,k:m),v(k:n,1));
A(k+1:n,k) = v(k+1:n,1);
end
function [v] = HOUSEHOLDER(a);
% input: a - vector with a6=0
% output: v - Householder vector of x
% such that v1 = 1
n = length(a);
v = a;
if (a(1) >= 0) beta = a(1) + norm(a);
else beta = a(1) - norm(a);
end
v(2:n) = 1/beta * v(2:n);
v(1) = 1;
|
|
|
| |
|
Danke für die Antwort, Adlerbob!
Ich habe inzwischen den Eindruck, dass die Matrix einfach sehr schlecht
konditioniert ist (Werte liegen extrem weit auseinander). Weil ich im
Prinzip nur den Eigenvektor zum kleinsten Eigenwert benötige, bin ich
auf folgendes Verfahren umgestiegen.
Bestimmung des Eigenvektors zum kleinsten Eigenwert:
[mm] x_{0} [/mm] sei ein beliebiger (vom 0-Vektor verschiedener) Start-Vektor
[mm] A_{} [/mm] ist eine reguläre nxn-Matrix
dann ist:
[mm] x_{1}=A^{-1}x_{0}
[/mm]
[mm] x_{i+1}=A^{-1}x_{i}
[/mm]
Zur Reduzierung von Rundungsfehlern werden alle Vektoren immer auf
Länge 1 normiert!
Bedingung für Iterationsende:
Wenn sich zwischen den Ergebnisvektoren der letzten beiden Iterationen
keine nennenswerte Änderung ergibt, kann abgebrochen werden. Keine
Veränderung heißt, die beiden letzten Vektoren sind parallel. Parallel
bedeutet im [mm] R^{2} [/mm] bzw. [mm] R^{3}, [/mm] dass der Winkel zwischen ihnen 0° ist - sein
Cosinus ist 1. Der Cosinus ergibt sich aus dem Skalarprodukt der (auf
Länge 1 normierten) Vektoren [mm] (). [/mm] Das gleiche Prinzip kann allge-
mein auf Vektoren im [mm] R^{n} [/mm] angewandt werden.
Das Verfahren lässt sich leicht zur Bestimmung des Eigenvektors zum
größten Eigenwert umformen:
[mm] x_{i+1} [/mm] = [mm] Ax_{i}
[/mm]
THEO
|
|
|
|
