A=[83.567	51.789;
51.789	51.344];
% Trovo il polinomio caratteristico
format bank
disp(poly(A))

% Autovalori ordinati in senso crescente
[Vini,Lambdaini]=eig(A);
lam=diag(Lambdaini);
[~,pos]=sort(lam,'descend');
% Autovettori e autovalori ordinati
V=Vini(:,pos);
Lambda=diag(lam(pos));

% Controllo che la traccia della matrice A sia uguale alla somma degli
% autovalori
assert(abs(trace(A)-sum(lam))<1e-12,"Errore di programmazione nella" + ...
    "verifica che la somma degli autovalori sia uguale alla traccia")

% Controllo che il determinante della matrice A sia uguale alla somma degli
% autovalori
assert(abs(det(A)-prod(lam))<1e-12,"Errore di programmazione nella" + ...
    "verifica che il prodotto degli autovalori sia uguale al determinante")
 
%% Verifica della scomposizione spettrale
Aresidua=A-V*Lambda*V';
assert(max(abs(Aresidua),[],'all')<1e-12,"Errore di programmazione nella" + ...
    "ricostruzione della matrice di partenza")


%% Verifo che se A viene ricostruita con la miglior matrice di rango 1
% la somma dei quadrati dei residui è il quadrato del secondo autovalore

Aresidua1=A-V(:,1)*Lambda(1,1)*V(:,1)';
% Somma dei quadrati dei residui tra A e la sua ricostruzione di rango 1
sumres2=sum((Aresidua1).^2,"all");
disp("Somma dei quadrati dei residui tra A e la sua ricostruzione di rango 1")
disp(sumres2)
disp("Secondo autovalore al quadrato")
disp(Lambda(2,2)^2)

assert(abs(sumres2-Lambda(2,2)^2)<1e-12,"Errore di programmazione nella" + ...
    "verifica della somma dei quadrati dei residui nella ricostruzione " + ...
    "della matrice di partenza con la migliore di rango 1")