Verfasst am: 01.11.2012, 13:28
Titel: Hilfe: Reordering Schur Form
Hallo,
ich habe ein Problem mit der Neuordnung der realen Schurform.
Und zwar benötige ich "[Q,R] = ordschur(Q,R,'lhp')"... das Programm was ich habe spuckt allerdings "[Q,R] = ordschur(Q,R,'udi')" aus. Leider finde ich nicht den Punkt, an dem ich es ändern muss - geschweige denn, dass ich weiß wie ich es ändern muss.
Meines Erachtens soll der BubbleSort verwendet werden, jedoch ist es so, dass das größte Element unten rechts in der Diagonalen stehen soll und das kleinste oben links.
%INPUT: Q...orthogonal, real % R...quasi-oberer Dreiecksmatrix, real - s.d. AQ = QR % z...Speicher in komplexer Ebene % b...bestimmt Länge der Ordnung
%
%
% b < 0, dann -b Blocks werden sortiert % b > 0, dann b oder b+1 Eigenwerte werden sortiert (abhängig von der % Größe der Blöcke % b = 0, komplette real Schurform wird sortiert
%
% lqr(A,B,Q,R) -> lösung überprüfen -> = newton %
r = find(abs(diag(R,-1)) > 100*eps); %Suche Einträge ungl. Null auf der Nebendiagonale
s = 1:size(R,1)+1; %Bilde Vektor mit der jeweils oberen linken Position von jedem Block
s(r+1) = [];
for k=1:length(s)-1;
sk = s(k);
if s(k+1) - sk == 2 %wenn block 2x2 -> normalize
[Q,R] = normalize(Q,R,sk:s(k+1)-1);
p(k) = R(sk,sk) + sqrt(R(sk+1,sk)*R(sk,sk+1)); % speichern der Eigenwerte else
p(k) = R(s(k),s(k)); %wenn block 1x1
end end
%richtige eigenwerte auswählen (nur negative in menge lambda)
for k = swaplist(p,s,z,b); %für k über alle nachbar-swaps
v = s(k):s(k+1)-1; %koordinaten der blocks
w = s(k+1):s(k+2)-1;
nrA = norm(R([v,w],[v,w]),inf); % berechne norm der matrix A der sylvestergleichung [Q,R] = swap(Q,R,v,w); %swap der blocks
s(k+1) = s(k)+s(k+2)-s(k+1); %update position
v = s(k):s(k+1)-1; %update koordinaten
w = s(k+1):s(k+2)-1;
iflength(v) ==2 [Q,R] = normalize(Q,R,v);
end iflength(w) == 2 [Q,R] = normalize(Q,R,w);
end
ap(k) = norm(R(w,v),inf)/(10*eps*nrA); %größe des blocks linksunten
end
R = R - tril(R,-2); %nullmatrix unterhalb einträge
for j = 2:length(s)-1; %erneute quasi-obere dreiecksmatrix
R(s(j), s(j)-1) = 0;
end
function[U,S] = normalize(U,S,v);
Q = rot(S(v,v)); %bestimme Givens-Rotation für die Standardisierung
S(:,v) = S(:,v)*Q; % links & rechts -> S, rechts -> U
S(v,:) = Q'*S(v,:); % nur Zeilen und Spalten mit Indizes im Vektor können betroffen sein
U(:,v) = U(:,v)*Q;
c = 1; %Start: Identität transformation
s = 0;
if X(1,1)~= X(2,2); %wenn nötig ändern zu...
tau = (X(1,2)+X(2,1))/(X(1,1)-X(2,2));
off = sqrt(tau^2+1);
v = [tau - off, tau +off];
[d,w] = min(abs(v));
c = 1/sqrt(1+v(w)^2);
s = v(w)*c;
end
Q = [c+s,s c];
function v = swaplist(p,s,z,b); %anders aufbauen, sodass ordschur(Q,R,'lhp') rauskommt
n = length(p);
k = 0;
v = [];
srtd = 0; %anzahl der sortierten eigenwerte
q = diff(s); %block größe
fini = 0;
while ~fini
k = k+1;
[dum,j] = select(p(k:n),z); %bestimme welcher block zur position k soll
p(k:n+1) = [p(j+k-1) p(k:n)]; %einfügen in k
p(j+k) = []; %löschen vom alten platz
q(k:n+1) = [q(j+k-1) q(k:n)]; %gleiche für block-größe
q(j+k) = [];
v = [v,j+k-2:-1:k]; %update liste der swaps für den block
srtd = srtd + q(k); %update anzahl der sortierten eigenwerte
fini = (k==n-1)|(k==-b)|(srtd == b)| ((srtd == b+1)&(b~=0));
end
function[val,pos] = select(p,z); % speicher zur oberen halbebene bewegen
y = real(z) +abs(imag(z))*i;
[val,pos] = min(abs(p-y)); %Finde nächsten block zum speicher
function[U,S] = swap(U,S,v,w);
[p,q] = size(S(v,w)); %p,q block größen
Ip = eye(p);
Iq = eye(q);
r = [];
for j = 1:q %vektorisiere für rechtsseitiges kronecker produkt
r = [r;S(v,w(j))]; %formuliere der sylvester gleichung
end
K = kron(Iq,S(v,v))-kron(S(w,w)',Ip); %kronecker produkt system matrix
[L,H,P,Q] = lu_complpiv(K); %LU-Zerlegung der matrix
gamma = min(abs(diag(H))); %scalierungs faktor um overlow vorzubeugen
sigp = 1:p*q;
for k = 1:p*q-1; %implementierung der permutation P der LU-Zerlegung PAQ = LU
sigp([k,P(k)]) = sigp([P(k),k]);
end
r = gamma*r(sigp); % skaliere und permutiere rechtseitig
x = (H\(L\r)); %löse die zwei dreieckssysteme
sigq = 1:p*q;
for k = 1:p*q-1;
sigq([k,Q(k)]) = sigq([Q(k),k]); %implementierung der permutation Q der LU-Zerlegung PAQ = LU end
x(sigq) = x; %permutiere die lösung
X = [];
for j = 1:q %entvektorisiere die lösung zurück in einen block oder beende die kronecker formulieren
X = [X,x((j-1)*p+1:j*p)];
end [Q,R] = qr([-X;gamma*Iq]); %Householder QR-Zerlegung von X
S(:,[v,w]) = S(:,[v,w])*Q;
S([v,w],:) = Q'*S([v,w],:);
U(:,[v,w]) = U(:,[v,w])*Q;
function[L,U,P,Q] = lu_complpiv(A); %LU-Zerlegung mit pivoting
P = [];
Q = [];
n = size(A,1);
for k = 1:n-1;
[a,r] = max(abs(A(k:n,k:n)));
[dummy,c] = max(abs(a));
cl = c+k-1;
rw = r(c)+k-1;
A([k,rw],:) = A([rw,k],:);
A(:,[k,cl]) = A(:,[cl,k]);
P(k) = rw;
Q(k) = cl;
if A(k,k) ~= 0;
rs = k+1:n;
A(rs,k) = A(rs,k)/A(k,k);
A(rs,rs) = A(rs,rs)-A(rs,k)*A(k,rs);
end end
U = tril(A')';
L = tril(A,-1)+ eye(n);
Ich denke, dass ich etwas in der Funktion swaplist ändern muss.
Es wäre super, wenn mir jemand helfen könnte - es muss auch nicht die Lösung sein, ein Tipp würde mir vielleicht auch schon weiter helfen.
Wenn jemand noch mehr Informationen braucht nur zu! :)
Ansonsten schon einmal DANKE! ;)
Einstellungen und Berechtigungen
Du kannst Beiträge in dieses Forum schreiben. Du kannst auf Beiträge in diesem Forum antworten. Du kannst deine Beiträge in diesem Forum nicht bearbeiten. Du kannst deine Beiträge in diesem Forum nicht löschen. Du kannst an Umfragen in diesem Forum nicht mitmachen. Du kannst Dateien in diesem Forum posten Du kannst Dateien in diesem Forum herunterladen
MATLAB, Simulink, Stateflow, Handle Graphics, Real-Time Workshop, SimBiology, SimHydraulics, SimEvents, and xPC TargetBox are registered trademarks and The MathWorks, the L-shaped membrane logo, and Embedded MATLAB are trademarks of The MathWorks, Inc.
RSS