power.m

% 
% Solution of Eigen-value problem using the Power Method
%

A = [2 -1 0 0;-1 2 -1 0;0 -1 2 -1; 0 0 -1 2] % coef matrix
y = [1; 0; 0; 0]                             % starting vector

fprintf('\nPower Method Demonstration\n\n');
highlow = input('Determine highest (1) or lowest (2) eigenvalue: ');

if highlow == 2  % lowest eigenvalue case: calc matrix inverse
  A = inv(A);
end;

lambda_old = 1;
error = 1;
count = 0;

while abs(error) > 0.01
  count = count + 1;

  y = A * y;          % new eigenvalue
  lambda = max(y);    % new lambda
  y = y/lambda;       % normalize

  error = (lambda - lambda_old)/lambda * 100;
  lambda_old = lambda;
  fprintf('# %i Lambda = %f Error = %f\n',count,lambda, error);
end;

if highlow == 1
  fprintf('\nHighest eigenvalue is %f\n',lambda);
else
  fprintf('\nLowest eigenvalue is %f\n',1/lambda);
end;