%% TMA671 Lösning av uppgift N.5.50
clear;
% Under storgruppsövningen som hölls onsdagen den 9e Maj 2018
% så demonstrerades uppgift 50 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 startgissningen
A = [6, 2, 1; 2, 3, 1; 1, 1, 1];% Matrisen
x = input('Startpunkt x-koordinat:');
y = input('Startpunkt y-koordinat:');
z = input('Startpunkt z-koordinat:');
x0 = [x; y; z]; % Startgissning
nuOfIter = input('Antal iterationer:'); % Antal iterationer
% Allokera minne till vår egenvektor
% Vi vill byta plats på rad 2 och 3
P = [0, 0, 1; 1, 0, 0; 0, 1, 0];
% Shift variabeln
shift = 2;


% Skriv till användaren
fprintf('\n\n=================================================================\n\n');
fprintf('\tUPPGIFT N.5.50: Approximera ett specifikt egenvärde\n\n');
fprintf('=================================================================\n\n');
fprintf('\tInvers iteration med skift\n');
fprintf('-----------------------------------------------------------------\n');

%for outerIndex = 1:length(nuOfIterVec)
% Definiera vår faktorisering
[L,U] = lu(P*(A-shift*eye(3)));
% Börja om från början
xk = x0;
% 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;
    % Uppdatera vår shift variabel
    shift = xk'*(A*xk);
    % Definiera vår faktorisering
    [L,U] = lu(P*(A-shift*eye(3)));
end
% Beräkna egenvärdet
lambda = (xk')*(A*xk);
%end

%xk

fprintf('-----------------------------------------------------------------\n');
% Jämför med eig
[V,D] = eig(A);
% Det är det andra egenvärdet vi är intresserade av
v = V(:,2);
lambdaFacit = D(2,2);

fprintf('Egenvärde:\n');
fprintf('\tInvers\t\t\tMatlab\n');
fprintf('\t%0.10f\t\t%0.10f\n',lambda,lambdaFacit);
format long
fprintf('Egenvektor:\n');
[xk,v]
fprintf('=================================================================\n\n');