clear();

//Fonction permettant de construire la matrice diagonale
function [Delta]=matDelta(k,n)
    
    
    // Initialisation de la matrice Delta
    Delta=zeros(n,n);
    
    //Remplissage de la première "case" de la matrice
    Delta(1)=1;

    //Remplissage du reste de la matrice
    for i=2:n
       Delta(i,i)=k^(-1/(n-1));
       n=n-1;
    end
    
endfunction

// Construction de la matrice H
function [H]=matH(n)
    
    // Initialisation de la matrice H
    H=zeros(n,n);
    
    // Remplissage de la diagonale de 2
    for i=1:n
        H(i,i)=2;
    end

    // Remplissage des deux autres diagonales de -1
    for i=1:n-1
        H(i,i+1)=-1;
        H(i+1,i)=-1;
    end
    
endfunction

// Construction de A
function [A]=matA(H,Delta)
    
    // Initialisation des matrices
    A=zeros(n,n);
    U=zeros(n,n);
    S=zeros(n,n);
    V=zeros(n,n);
    
    // Récupération de U et V
    [U,S,V] = svd(H);
    
    //Calcul de A
    A=U*Delta*V';
    
endfunction

// Construction de Eps
function [Eps]=matEps(n,val)
    
    // Initialisation de la matrice Eps
    Eps=zeros(n,1);
    
    // Remplissage de la matrice de la valeur souhaité (ici 1)
    for i=1:n
        Eps(i,1)=val;
    end
    
endfunction

// Fonction pour la méthode du Jacobi
function [x]=Jacobi(k,n,D,L,U,b)

    if k==0 then
        x=zeros(n,1);
    elseif (k<=1 & k > 0)
        x=(D^-1)*b;
    else
        x= -(D^-1)*(L+U)*Jacobi(k-1,n,D,L,U,b) + (D^-1)*b;
    end

endfunction




// Utilisation des fonctions avec données du sujet
    // Entrée des données selon le calcul souhaité
        kappa=10; n=10;

    // Génération des matrices
        Delta = matDelta(kappa,n);
        H = matH(n);
        A = matA(H,Delta);
        Eps = matEps(n,1);
        b = A*Eps;
        disp(A), disp(b);
    // Calcul de la solution de Ax=b
        [x] = linsolve(A,-b);
        //disp("x="), disp(x);
        
    // Calcul des normes pour l'évaluation de l'erreur
        norme_x = norm(x);
        norme_Eps = norm(Eps);
        //disp("norme_x="), disp(norme_x); disp("norme_Eps="), disp(norme_Eps);
        
    // Récupération des différentes partie de A
        d = diag(A), D = diag(d);
        L = tril(A) - D;
        U = triu(A) - D;
        //disp(A); disp(D); disp(U); disp(L);
        
    // Utilisation de la méthode du Jacobi
        k=110;    // Fixation de k pour un seuil < 10^-8 (donc déterminé après test)
        xJac=Jacobi(k,n,D,L,U,b);
        disp("xJac="),disp(xJac);
        
        // Calcul de l'erreur résiduelle        
            r = (norm(A*xJac -b))/norm(b);
            e = norme_Eps - norme_x;
            disp("L''erreur résiduelle est égale à"), disp(r);
            disp("L''erreur par la méthode directe est"), disp(e);
        
        
        
        
        
        
        
        
        
        
        