%% Demonstrate Eigenvalue Calculations via Power Method
%
clear;

K=[12 -4 0; -4 6 -2; 0 -2 2];
M=[2 0 0; 0 2 0; 0 0 2];
A=K*inv(M);

%Let Matlab Calculate
display('Calculate eigenvectors (V) and eigenvalues (D) ...');
[V,D]=eig(A)
display('Calculate |A-lI| ... ');
det(A-D(1,1)*eye(3))

%Set iteration parameters
toler=0.0001;
maxiter=100;
n=size(A,1);

%Calculate by Power Method Iteration
display('Power Method Iteration');
xk=ones(n,1);
lambda0=1;
lambdak=0;
k=0;
while ((max(abs((lambdak-lambda0)/lambda0))>toler) & (k<maxiter));
    lambda0=lambdak;
    k=k+1;
    x0=xk;
    xk=A*x0;
    lambdak=xk(1)./x0(1);
    display(['k=',num2str(k,'%3.3i'),' | ',num2str(xk','  %1.16e'),' | ',num2str(lambdak,'  %1.16e')]);
end;

%Calculate by Scaled Power Method Iteration (Find Extreme, Largest)
display('Scaled Power Method Iteration');
xk=ones(n,1);
lambda0=1;
lambdak=0;
k=0;
while ((max(abs((lambdak-lambda0)/lambda0))>toler) & (k<maxiter));
    lambda0=lambdak;
    k=k+1;
    x0=xk;
    xk=A*xk;
    lambdak=xk(1)./x0(1);
    xk=xk./max(xk);
    display(['k=',num2str(k,'%3.3i'),' | ',num2str(xk','  %1.16e'),' | ',num2str(lambdak,'  %1.16e')]);
end;

%Calculate by Shifted Scaled Power Method Iteration (Find Opposite Extreme)
display('Shifted Scaled Power Method Iteration');
shift=lambdak;
xk=ones(n,1);
lambda0=1;
lambdak=0;
k=0;
while ((max(abs((lambdak-lambda0)/lambda0))>toler) & (k<maxiter));
    lambda0=lambdak;
    k=k+1;
    x0=xk;
    xk=(A-shift*eye(n))*xk;
    lambdak=xk(1)./x0(1);
    xk=xk./max(xk);
    display(['k=',num2str(k,'%3.3i'),' | ',num2str(xk','  %1.16e'),' | ',num2str(shift+lambdak,'  %1.16e')]);
end;
