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
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