using data from while loop and resetting initial conditions

I am writing an algorithm for the power method, which is an iterative method for solving for all eigenvalues and eigenvectors of a matrix. Finding the first set of data, the dominant eigenvalue/vector is fairly straightforward.

What we do next is where i get confused. we are supposed to take the data and apply the following formula:. We then use Bv in the next iteration. I do not know how to take this value and input it into the new iteration.

The code:

clc;
%Assignement Conditions
    FA = @(k_1, k_2, k_3,m_1, m_2) [-(k_2+k_1)/m_1 , k_2/m_1 ; k_2/m_2 -(k_3+k_2)/m_2];
    A = FA(1,1,1,2,1);

%Stopping Criteria
    e_s = .00005;
    disp(['<strong> Stopping Criteria: </strong>' strcat(num2str(e_s),'%') ])
    
%Vector to display all Eigenvalues
    Eigenvalue_summary = zeros(1,height(A));

%Dominant Eigenvalue 

    %Random Guess for Dominant Eigenvector 
    V_approximation = ones(height(A),1,'double');

    %Iterative method for Dominant Eigenvalue/Eigenvector
    Iteration = 0;
    while (1)
    disp('====================')
    Iteration = Iteration + 1;
    disp(['<strong>Iteration: </strong>' num2str(Iteration) ])
    
    %Displaying Eigenvector to be used
    disp(['<strong>Initial Guess: </strong>'])
    for i = 1:length(V_approximation)    
    fprintf('%f\n',V_approximation(i))
    end
    

    %Sign of previous V (before iteration)
    sign_old = sign(V_approximation);

    %Approximations
    Eigenvalue = norm(A*V_approximation);
    Eigenvector = (A*V_approximation)/Eigenvalue;

    %Sign of new V
    sign_new = sign(Eigenvector);



    %Reseting Eigenvalue based on sign of current and prev. V
    if sign_new ~= sign_old
    Eigenvalue = -Eigenvalue;
    end


    %Recording Approximations 
    for count = Iteration
    C{count} = table(Iteration,Eigenvalue);
    D{count} = table([Eigenvector]);
    end

    Eigen_Val = vertcat(C{:});
    Eigen_Vec = vertcat(D{:});

    if Iteration > 1
    error = string(abs((table2array(Eigen_Val(Iteration,2)) - table2array(Eigen_Val(Iteration-1,2)))/table2array(Eigen_Val(Iteration,2)))*100);
    else
    error = string('----');
    end



    for count = Iteration
    E{count} = table(Iteration,Eigenvalue,error);
    end
    Error = vertcat(E{:});

    disp(['<strong>Eigenvalue: </strong>' num2str(Eigenvalue)])
    
    
    disp(['<strong>Eigenvector: </strong>'])
    for i = 1:length(Eigenvector)    
    fprintf('%f\n',Eigenvector(i))
    end

    disp(['<strong>Relative Error: </strong>' num2str(error) ])
    
    
    %Resetting Eigenvector Approximation for next iteration
    V_approximation = Eigenvector;

    if str2num(error) < e_s,break,end

    end

    %Tables of Iteration Data
    Eigen_Vec
    Error
    
    %Final Results
    Dominant_Eigenvector = Eigenvector
    Dominant_Eigenvalue = Eigenvalue
    Eigenvalue_summary(1) = Dominant_Eigenvalue;



    B1 = A - Eigenvalue_summary(1)*eye(height(A))

Note the final line of code, where B is determined

Thank you

Answer 1

[Chetan]

Hi @Michael Jacobson,

I understand that you're implementing the power method in MATLAB to find the dominant eigenvalue and eigenvector of a matrix. And wish to modify the matrix for subsequent iterations to find more eigenvalues and eigenvectors.

You can modify the matrix by subtracting the outer product of the found eigenvector, scaled by the found eigenvalue, from the original matrix. This is known as deflation.

Here's how you can adjust your code to include this step:

clc;
% Assignment Conditions
FA = @(k_1, k_2, k_3, m_1, m_2) [-(k_2+k_1)/m_1 , k_2/m_1; k_2/m_2, -(k_3+k_2)/m_2];
A = FA(1, 1, 1, 2, 1);

% Stopping Criteria
e_s = 0.00005;
disp(['<strong> Stopping Criteria: </strong>', num2str(e_s), '%'])

% Vector for Eigenvalues
Eigenvalue_summary = zeros(1, height(A));

% Dominant Eigenvalue
V_approximation = ones(height(A), 1, 'double'); % Initial guess

Iteration = 0;
while true
    Iteration = Iteration + 1;
    disp(['<strong>Iteration: </strong>', num2str(Iteration)])
    
    Eigenvalue = norm(A*V_approximation);
    Eigenvector = (A*V_approximation) / Eigenvalue;
    
    % Update for next iteration
    V_approximation = Eigenvector;
    
    % Check convergence
    if Iteration > 1 && abs(Eigenvalue - Eigenvalue_summary(Iteration-1)) < e_s
        break;
    end
    Eigenvalue_summary(Iteration) = Eigenvalue;
    
    % Modify A for next eigenvalue (Deflation)
    A = A - Eigenvalue * (Eigenvector * Eigenvector');
end

% Final Results
Dominant_Eigenvector = Eigenvector;
Dominant_Eigenvalue = Eigenvalue;
disp(['Dominant Eigenvalue: ', num2str(Dominant_Eigenvalue)]);
disp('Dominant Eigenvector: ');
disp(Dominant_Eigenvector);

This code iteratively finds the dominant eigenvalue and eigenvector, then deflates the matrix `A` for subsequent iterations.

Refer to the following MathWorks Documentation for more details:

• https://www.mathworks.com/help/matlab/ref/eig.html

Thanks

Chetan