%% TMA671 Lösning av uppgift N.5.27
% Under storgruppsövningen som hölls onsdagen den 9e Maj 2018
% så demonstrerades uppgift 27 i Kapitel 5 i Numerisk Analys
% boken. Uppgiften handlar om att approximera egenvärdena till 
% en matris med hjälp av potensmetoden och inversiteration
% Programmet är skrivet av Johannes Borgqvist

%% Definiera matrisen och startuppgiften
A = [2, 1; 2, 3];% Matrisen
x0 = [1; 1]; % Startgissning
nuOfIter = input('Ange antal iterationer:'); % Antal iterationer
% Allokera minne till vår egenvektor
xk = x0;
% Skriv ut ett meddelande till användaren
fprintf('=================================================================\n\n');
fprintf('\tUPPGIFT N.5.27: Approximera egenvärden\n\n');
fprintf('=================================================================\n\n');

%% Delupgift (a): beräkna egenvärden med potensmetoden
fprintf('-----------------------------------------------------------------\n');
fprintf('\tDeluppgift (a): Potensmetoden\n');
fprintf('-----------------------------------------------------------------\n');
% Loopa över antalet iterationer
for index = 1:nuOfIter
    % Skriv ut iterationen
    fprintf('Iteration %d\n',index);
    % Börja iterera
    yk = A*xk; 
    xk = (1/norm(yk))*yk 
end
% Beräkna egenvärdet
lambdaPotens = (xk')*(A*xk);
fprintf('\nPotens egenvärde:\t\tlamdbda\t=\t%0.5f\n\n',lambdaPotens);


%% Delupgift (b): beräkna egenvärden med invers iteration
fprintf('-----------------------------------------------------------------\n');
fprintf('\tDeluppgift (b): Invers iteration\n');
fprintf('-----------------------------------------------------------------\n');
% Börja om från början
xk = x0;
% Definiera permuteringsmatris 
P = [0, 1; 1, 0];
% Vår L matris 
L = [1, 0; 1, 1,];
% Vår U matris
U = [2, 3; 0, -2];
% Okaj, nu loopar vi över antal iterationer
for index = 1:nuOfIter
    % Skriv ut iterationen
    fprintf('Iteration %d\n',index);
    % Börja iterera: lös för z
    zk = L\(P*xk); 
    % Lös för y
    yk = U\zk;
    % Definiera xk
    xk = (1/norm(yk))*yk 
end
% Beräkna egenvärdet
lambdaInvers = (xk')*(A*xk);
fprintf('\nInvers egenvärde:\t\tlamdbda\t=\t%0.5f\n\n',lambdaInvers);
fprintf('=================================================================\n\n');