Home > K25104 > LinearSystems > inverse_iteration.m

inverse_iteration

PURPOSE ^

Implements inverse iteration to find an eigenvalue of A

SYNOPSIS ^

function [ v,c,k ] = inverse_iteration( A,s,tol )

DESCRIPTION ^

 Implements inverse iteration to find an eigenvalue of A 
   and its corresponding eigenvector
 Input arguments:
   A, square matrix
   s, estimate of eigenvalue
   tol, tolerance
 Output arguments:
   v, eigenvector of A
   c, corresponding eigenvalue
   k, number of iterations

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

0001 function [ v,c,k ] = inverse_iteration( A,s,tol )
0002 % Implements inverse iteration to find an eigenvalue of A
0003 %   and its corresponding eigenvector
0004 % Input arguments:
0005 %   A, square matrix
0006 %   s, estimate of eigenvalue
0007 %   tol, tolerance
0008 % Output arguments:
0009 %   v, eigenvector of A
0010 %   c, corresponding eigenvalue
0011 %   k, number of iterations
0012 
0013 [n,m]=size(A); % finding the size of A
0014 if n~= m;
0015     error('input is not a square matrix');
0016 elseif tol<=0;
0017     error('tolerance should be positive');
0018 end
0019 
0020 k=0;    % initializing k
0021 
0022 while k>=0; 
0023     [L,U,P]=LU_pivot(A-s*eye(n));   % calculate the LU factorization
0024                                     % with pivoting of (A-sI)
0025     d=diag(U);      
0026     if min(abs(d))<=1e-12;  % if any diagonal elements of U are close to zero
0027         c=s;                % accept s as the eigenvalue
0028         v=U;                % set output v to U
0029         disp(['eigenvalue found; v lies in the null space of '...
0030             'given matrix']);
0031         return;             % return function
0032     else
0033         if k==0;
0034             [~,i]=min(abs(d));      % i is the index of the smallest entry
0035                                     % along the diagonal of U
0036             ei=zeros(n,1); ei(i)=1; % initialize ith standard basis vector
0037             x0=L*ei;                % initialize x0
0038             x1=Backward(U,ei);     % calculate x1
0039             x1=P*x1; x0=P*x0;   % apply pivoting
0040         else
0041             y=Forward(L,x0);      % calculate x1 from (A-sI)x1=x0
0042             x1=Backward(U,y);  
0043             x1=P*x1;            % apply pivoting
0044         end
0045         
0046         mu=x1'*(A-s*eye(n))*x1/(x1'*x1); % calculate mu which minimizes
0047                                          % |x0-mu*x1|
0048         z=x0-mu*x1;            % if |z|<tol|x1| we terminate procedure
0049         if sqrt(z'*z)<tol*sqrt(x1'*x1); 
0050             c=mu+s;            % accept mu+s as eigenvalue
0051             v=x1/sqrt(x1'*x1);  % normalize x1 and accept as eigenvector
0052             return;     % return function
0053         else
0054             x0=x1/sqrt(x1'*x1); % normalize x
0055             s=s+mu;     % adjust s
0056             k=k+1;      % increment k
0057         end
0058     end
0059 end
0060 end
0061             
0062             
0063

Generated on Mon 18-Jan-2016 10:25:49 by m2html © 2005