Lowest Positive Eigenvalue using the inverse iteration method

eigenvalues-eigenvectors

I am trying to find the lowest positive eigenvalue of the generalized eigenproblem of buckling using the inverse iteration method. I have successfully implemented the inverse iteration to find the absolute lowest eigenvalue. However, i do not know how to perform the shift on the matrices to find a positive eigenvalue.

This is my code so far in MatLab used for testing:

%Generalized EigenProblem: Ke.x + λ.Kg.x = 0
Ke=[1,2,-2;2,4,2;-2,2,4];
Kg=[0,0,3;6,7,2;4,4,5];

%Matlab Result
d=eigs(Ke,-Kg);
fprintf('Matlab EigenValues:  %f, %f, %f \n',d(1),d(2),d(3))

%Initial Run to Find the Smallest Absolute Value 
Err = 10;
EigenValueEstimate =1;
EigenVectorEstimate =[1;1;1];
y1 =[0;0;0];
y2 = [0;0;0];
TotalIterations = 0;
while Err > 10^-5 && TotalIterations <= 100
    TotalIterations = TotalIterations + 1;
    %Solve
    y1=-Kg*EigenVectorEstimate;
    EigenVectorEstimate=Ke\y1;
    %Norm
    y2=-Kg*EigenVectorEstimate;
    normV  = EigenVectorEstimate'* y2;
    EigenVectorEstimate= 1 / sqrt(abs(normV)).* EigenVectorEstimate;
    %EigenValue
    y2=Ke*EigenVectorEstimate;
    NewEigenValueEstimate = y2'* EigenVectorEstimate*(normV / abs(normV));
    Err = abs((NewEigenValueEstimate - EigenValueEstimate) / EigenValueEstimate);
   EigenValueEstimate = NewEigenValueEstimate;
end
fprintf('Lowest Absolute EigenValue is %f\n',EigenValueEstimate)

%Shift
Shift=-EigenValueEstimate*2;
Ke=Ke-eye(3)*(Shift);
Kg=Kg+eye(3)*(1/Shift);
%d=eigs(Ke,-Kg)

%ReRun
Err = 10;
EigenValueEstimate2 =1;
EigenVectorEstimate2 =[1;1;1];
y1 =[0;0;0];
y2 = [0;0;0];
TotalIterations = 0;
while Err > 10^-5 && TotalIterations <= 100
    TotalIterations = TotalIterations + 1;
    %Solve
    y1=-Kg*EigenVectorEstimate2;
    EigenVectorEstimate2=Ke\y1;
    %Norm
    y2=-Kg*EigenVectorEstimate2;
    normV  = EigenVectorEstimate2'* y2;
    EigenVectorEstimate2= 1 / sqrt(abs(normV)).* EigenVectorEstimate2;
    %EigenValue
    y2=Ke*EigenVectorEstimate2;
    NewEigenValueEstimate2 = y2'* EigenVectorEstimate2*(normV / abs(normV));
    Err = abs((NewEigenValueEstimate2 - EigenValueEstimate2) / EigenValueEstimate2);
    EigenValueEstimate2 = NewEigenValueEstimate2;
end
fprintf('Lowest Positive Absolute EigenValue is %f\n\n',EigenValueEstimate2+Shift)

Best Answer

The original EigenProblem is:

Ke.x + λ.Kg.x = 0

The same EigenProblem can be shifted by μ as follows:

(Ke - Kg.μI).x + λ.Kg.x = 0

The method to find the lowest positive eigenvalue is as follows:

  • Run Inverse Iteration Method to arrive at the lowest absolute eigenvalue
  • If the calculated eigenvalue is positive, accept the result
  • If it is negative, perform a shift by the first eigenvalue and run the inverse iteration again
  • Keep performing last step until a positive eigenvalue is arrived at

Below is the updated MatLab code where i am shifting by twice the first eigenvalue.

%Generalized EigenProblem: Ke.x + λ.Kg.x = 0
Ke=[1,2,-2;2,4,2;-2,2,4];
Kg=[0,0,3;6,7,2;4,4,5];

%Matlab Result
[d,p]=eigs(Ke,-Kg);
fprintf('Matlab EigenValues:  %f, %f, %f \n',p(1),p(5),p(9))

%Initial Run to Find the Smallest Absolute Value 
Err = 10;
EigenValueEstimate =1;
EigenVectorEstimate =[1;1;1];
y1 =[0;0;0];
y2 = [0;0;0];
TotalIterations = 0;
while Err > 10^-7 && TotalIterations <= 100
    TotalIterations = TotalIterations + 1;
    %Solve
    y1=-Kg*EigenVectorEstimate;
    EigenVectorEstimate=Ke\y1;
    %Norm
    y2=-Kg*EigenVectorEstimate;
    normV  = EigenVectorEstimate'* y2;
    EigenVectorEstimate=(1 / sqrt(abs(normV))).* EigenVectorEstimate;
    %EigenValue
    y2=Ke*EigenVectorEstimate;
    NewEigenValueEstimate = y2'* EigenVectorEstimate*(normV / abs(normV));
    Err = abs((NewEigenValueEstimate - EigenValueEstimate) / EigenValueEstimate);
    EigenValueEstimate = NewEigenValueEstimate;
end
fprintf('Lowest Absolute EigenValue is %f\n',EigenValueEstimate)

%Shift
Shift=-2*EigenValueEstimate;
Ke=Ke+Kg*(eye(3).*Shift);

%ReRun
Err = 10;
EigenValueEstimate2 =1;
EigenVectorEstimate2 =[1;1;1];
TotalIterations = 0;
while Err > 10^-7 && TotalIterations <= 100
    TotalIterations = TotalIterations + 1;
    %Solve
    y1=-Kg*EigenVectorEstimate2;
    EigenVectorEstimate2=Ke\y1;
    %Norm
    y2=-Kg*EigenVectorEstimate2;
    normV  = EigenVectorEstimate2'* y2;
    EigenVectorEstimate2= 1 / sqrt(abs(normV)).* EigenVectorEstimate2;
    %EigenValue
    y2=Ke*EigenVectorEstimate2;
    NewEigenValueEstimate2 = y2'* EigenVectorEstimate2*(normV / abs(normV));
    Err = abs((NewEigenValueEstimate2 - EigenValueEstimate2) / EigenValueEstimate2);
    EigenValueEstimate2 = NewEigenValueEstimate2;
end
fprintf('Lowest Positive Absolute EigenValue is %f\n\n',EigenValueEstimate2+Shift)
Related Question