Home > K25104 > LinearSystems > twostage_power.m

twostage_power

PURPOSE ^

Implements 2-stage power method to find a complex conjugate pair

SYNOPSIS ^

function [ v1,c1,v2,c2,r1,r2,k ] = twostage_power( A,x0,tol )

DESCRIPTION ^

 Implements 2-stage power method to find a complex conjugate pair 
   of eigenvalues for A, and their corresponding eigenvectors.
 Input arguments:
   A, square matrix
   x0, starting vector
   tol, tolerance
 Output arguments:
   v1,v2, eigenvectors of A
   c1,c2, corresponding eigenvalues
   r1,r2, real vectors which span the eigenspace 
   k, number of iterations

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

0001 function [ v1,c1,v2,c2,r1,r2,k ] = twostage_power( A,x0,tol )
0002 % Implements 2-stage power method to find a complex conjugate pair
0003 %   of eigenvalues for A, and their corresponding eigenvectors.
0004 % Input arguments:
0005 %   A, square matrix
0006 %   x0, starting vector
0007 %   tol, tolerance
0008 % Output arguments:
0009 %   v1,v2, eigenvectors of A
0010 %   c1,c2, corresponding eigenvalues
0011 %   r1,r2, real vectors which span the eigenspace
0012 %   k, number of iterations
0013 
0014 [n,m]=size(A); % finding the size of A
0015 if n~= m;
0016     error('input A is not a square matrix');
0017 elseif tol<=0;
0018     error('tolerance should be positive');
0019 end
0020 
0021 x1=A*x0;    % initialize x1
0022 k=0;        % initialize k to zero
0023 while k>=0;     % this will always hold, another statement within the loop
0024                 % determines convergence and returns the function
0025                 
0026     x2=A*x1;    % calculate x2
0027     
0028     % calculate a and b which minimize the norm of x2+a*x1+b*x0
0029     denom=((x0'*x0)*(x1'*x1)-(x1'*x0)^2);   
0030     a=(x1'*x0)*(x2'*x0)-(x0'*x0)*(x2'*x1);  
0031     a=a/denom;
0032     b=(x1'*x0)*(x2'*x1)-(x1'*x1)*(x2'*x0);
0033     b=b/denom;
0034     
0035     y=x2+a*x1+b*x0;  
0036     if norm(y)<tol;         % if |y|<tol we terminate procedure
0037         p=[1,a,b];          % then we solve the polynomial with these
0038         q=roots(p);         % coefficients
0039         
0040         c1=q(1); c2=q(2);   % take roots the polynomial to be the required
0041                             % eigenvalues
0042                             
0043         v1=x1-c2*x0;        % generate the corresponding eigenvectors
0044         v2=x1-c1*x0;
0045         
0046         r1=real(v1);        % let r1 and r2 be the real and imaginary parts
0047         r2=imag(v1);        % of the eigenvectors, then they are real
0048                             % vectors which span the same eigenspace as
0049                             % v1 and v2.
0050         return;             % return the function
0051     else
0052         x0=x1/norm(x1);     % adjust x0 and normalize
0053         x1=x2/norm(x1);     % adjust x1
0054         k=k+1;              % increment k
0055     end
0056 end
0057 end
0058 
0059 
0060     
0061

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