function lambda = my_eig(A,type) if nargin ==1; type=0; end [m,m] = size(A); lambda = zeros(m,1); for k = 1:m if (type==0) [lambda(k),A] = QR_rayleigh(A); elseif (type==1) [lambda(k),A] = QR_wilkinson(A); else disp('Unzulässiger Wert für type') end end return; % ------------------------------------------------- function [lambda,A_defl] = QR_rayleigh(A) max_steps = 10; I = eye(size(A)); norm_A = norm(A,inf); for steps = 0:max_steps mu = rayleighQuotient(A); [Q,R] = qr(A-mu*I); A = R*Q+mu*I; if (norm(A(end,1:end-1)) < eps*norm_A) break; end end lambda = A(end,end); A_defl = A(1:end-1,1:end-1); % ------------------------------------------------- function mu = rayleighQuotient(A) mu = A(end,end); % ------------------------------------------------- function [lambda,A_defl] = QR_wilkinson(A) max_steps = 10; [m,m] = size(A); I = eye(size(A)); norm_A = norm(A,inf); if (m==1) % 1d-Fall ist trivial A_defl=[]; lambda=A; return; end for steps = 0:max_steps mu = wilkinson(A); [Q,R] = qr(A-mu*I); A = R*Q+mu*I; if (norm(A(end,1:end-1)) < eps*norm_A) break; end end lambda = A(end,end); A_defl = A(1:end-1,1:end-1); % ------------------------------------------------- function mu = wilkinson(A) a = A(end,end); x=eig(A(end-1:end,end-1:end)); % nimm EW der minimalen Abstand zu a besitzt if (abs(x(1)-a)