addpath("./functions");
pkg load symbolic

function [K, K_built] = cond_inf(alpha)

A = [4 0 1 0; 0 alpha-1 0 0; 1 0 1 0; 0 0 0 2];
[L, U, err] = gauss_simple(A);

% Per la risoluzione si sfrutta l'equazione matriciale A * A_inv = I,
% in particolare la forma A * A_inv(:, j) = I(:,j)

A_inv = zeros(4,4);
for j = 1:4
    % aggiornamento del termine noto per la j-esima colonna
    b = zeros(4,1);
    b(j) = 1;

    % in base al valore del termine noto si ottengono i valori
    % della j-esima colonna di A_inv.
    [x, err] = lusolve(L, U, [], [], b);
    A_inv(:,j) = x;
end

% calcolo indice di condizionamento
K = norm(A, inf) * norm(A_inv, inf);
K_built = cond(A, inf);

end

alpha = 1.2 : 0.05 : 3.2;
len_alpha = length(alpha);
K = zeros(len_alpha, 1);
K_built = K;

for i = 1:len_alpha
    [x, y] = cond_inf(alpha(i));
    K(i) = x;
    K_built(i) = y;
end

plot( alpha, K, "ro--", ...
      alpha, K_built, "b--", ...
      1.2, 5:30);
legend("cond LU", "cond built-in")
ylabel('K_{inf}(A)')
xlabel('\alpha');
box off

syms alpha
A = [[4 0 1 0]; 0 alpha-1 0 0; [1 0 1 0]; [0 0 0 2]];
disp('Determinante di A: ')
det(A)
disp('Valori che annullano il determinante: ')
solve(det(A))

A = double(subs(A, alpha, 2));
syms x
p_char = det(A - x*eye(4))
solve(det(A - x*eye(4)))
x3 = (-sqrt(13) + 5) / 2
x4 = (sqrt(13) + 5) / 2

der_p_char = diff(p_char)

p_char = @(x) x.^4 - 8 * x.^3 + 20 * x.^2 - 19 * x + 6;
der_p_char = @(x) 4 * x.^3 - 24 * x.^2 + 40 * x - 19;

figure
t = linspace(0, 5);
plot( t, p_char(t), ...
      t, der_p_char(t), ...
      1, p_char(1), "ro", ... % zeri di p_char
      2, p_char(2), "ro", ...
      x3, p_char(x3), "ro", ...
      x4, p_char(x4), "ro");
legend("polinomio caratteristico", "derivata");
box off

format long
[x3_newton, appr_x3, nit_x3] = newton(p_char, der_p_char, ...
    0, 10^-12, 10^-12, 100);
[x4_newton, appr_x4, nit_x4] = newton(p_char, der_p_char, ...
    5, 10^-12, 10^-12, 100);
x3_newton
x3
x4_newton
x4

fprintf("Numero di iterazioni richieste per approssimare x3: %.i\n", ...
    nit_x3);
fprintf("Numero di iterazioni richieste per approssimare x4: %.i\n", ...
    nit_x4);

n = (10 : 1 : 50);
N = length(n);
err_rel_chol = zeros(N, 1);
err_rel_gauss = zeros(N, 1);

for k = 1:N
    B = (1 / n(k)) * hilb(n(k)) + diag(1:n(k));
    A = B' * B;
    
    x_esatta = (1:n(k))';
    b = A * x_esatta;
    
    [R, p] = chol(A, 'lower');
    %x_chol = lusolve(R, R', [], [], b);
    y = lsolve(R, b);
    x_chol = usolve(R', y);
    
    [L, U, err] = gauss_simple(A);
    x_gauss = lusolve(L, U, [], [], b);
    
    err_rel_chol(k) = norm(x_esatta - x_chol) / ...
        norm(x_esatta);
    err_rel_gauss(k) = norm(x_esatta - x_gauss) / ...
        norm(x_esatta);
end

figure
plot( n, err_rel_chol, "-*", ...
      n, err_rel_gauss, "o-");
legend("err rel cholesky", "err rel gauss");
box off