///Exercice 1 
//question 1: fonctions descente et remontee
disp("exercice 1 question 1")
function x=remontee(A,b)
    [n,m]=size(A);
    x=zeros(b);
    x(n)=b(n)/A(n,n);
    for i=n-1:-1:1
        x(i)=b(i)
        for k=i+1:n
            x(i)=x(i)-A(i,k)*x(k);
        end
        x(i)=x(i)/A(i,i)
    end
endfunction

//test
A1=[1 2 1;0 1 1;0 0 2];
b1=[1;2;3]
x1=remontee(A1,b1)
 //on vérifie qu'on a bien obtenu le bon résultat
if x1==A1\b1, then disp("remontee fonctionne sur cet exmple") end

function x=descente(A,b)
    [n,m]=size(A);
    x=zeros(b);
    x(1)=b(1)/A(1,1);
    for i=2:n
        x(i)=b(i)
        for k=1:i-1
            x(i)=x(i)-A(i,k)*x(k);
        end
        x(i)=x(i)/A(i,i)
    end
endfunction
//test 

A2=[1 0 0;2 2 0;1 0 2];
b2=[1;2;3]
x2=descente(A2,b2)
if x2==A2\b2, then disp("descente fonctionne sur cet exemple") end

//Question 2
disp("exercice 1 question 2")
A=[2 1 1 0;4 3 3 1;8 7 9 5;6 7 9 8]
[n,m]=size(A)
mineurs=[]
test=[]
for i=1:n
    mineurs=[mineurs; det(A(1:i,1:i))]
    if det(A(1:i,1:i))==0,
        then test=[test; 0]
        else test=[test; 1] 
    end    
end
disp("les mineurs de A valent :")
disp(mineurs)
if (test==1) then
    disp("A admet une decompostion LU")
    else disp("A n admet pas de decomposition LU")
end


//Question 3 
disp("exercice 1 question 3")
[L,U]=lu(A)// Attention L n'est pas triangulaire inf !!!
disp("matrice L obtenue par Scilab avec [L,U]=lu(A)")
disp(L)
[L,U,E]=lu(A) // E matrice de permutation tq EA=LU
disp("matrice L obtenue avec [L,U,E]=lu(A)")
disp(L)
b=ones(4,1)
c=E*b
Z=descente(L,c)
X=remontee(U,Z)

//Exercice 2
disp("Exercice 2")
function L=cholesky(A)
    [n,m]=size(A);
    L=zeros(A);
    L(1,1)=sqrt(A(1,1))
    L(:,1)=A(:,1)/L(1,1);
    for j=2:n
        L(j,j)=sqrt(A(j,j)-sum(L(j,1:j-1).^2))
        for i=j+1:n
            L(i,j)=(A(i,j)-sum(L(i,1:j-1).*L(j,1:j-1)))/L(j,j);
        end
    end
endfunction

A=[1 1 -1 1;1 3 -1 5; -1 -1 2 -2; 1 5 -2 13];//A est bien symetrique
if (A'==A) then
    disp("A est symétrique")
    else disp("A n est pas symétrique")
end
spec(A) //on verifie que les valeurs propres sont >0
if (spec(A)>0) then
    disp("A est définie positive")
    else disp("A n est pas definie positive")
end
L=cholesky(A) //on verifie que L est bien triang inf et que LL'=A
disp("Matrice L calculée avec les formules du cours")
disp(L)
Q=chol(A) // on obtient la transposée de L !
disp("matrice Q obtenue avec Q=chol(A)")
disp(Q)



//Exercice 3 
disp("exercice 3")
//question 1
disp("question 1")
A=[1 2 -2; 1 1 1; 2 2 1];
disp("A=")
disp(A)
//Matrice de Jacobi 
MJ=diag(diag(A)); NJ=MJ-A;
rhoJ=max(abs(spec(inv(MJ)*NJ)))//rayon spectral de J
disp("rayon spectral de J")
disp(rhoJ)

[n,m]=size(A)
MG=zeros(A)
MG(1,1)=A(1,1)
for i=2:n
    MG(i,1:i-1)=-A(i,1:i-1);
    MG(i,i)=A(i,i);
end
NG=MG-A;
rhoG=max(abs(spec(inv(MG)*NG)))
disp("rayon spectral de GS")
disp(rhoG)
//Jacobi CV, GS DV

//question 2 
disp("question 1")
A=[2 -1 1; 2 2 2; -1 -1 2];
disp("A=")
disp(A)
//Matrice de Jacobi 
MJ=diag(diag(A)); NJ=MJ-A;
rhoJ=max(abs(spec(inv(MJ)*NJ)))//rayon spectral de J
disp("rayon spectral de J")
disp(rhoJ)

[n,m]=size(A)
MG=zeros(A)
MG(1,1)=A(1,1)
for i=2:n
    MG(i,1:i-1)=-A(i,1:i-1);
    MG(i,i)=A(i,i);
end
NG=MG-A;
rhoG=max(abs(spec(inv(MG)*NG)))
disp("rayon spectral de GS")
disp(rhoG)
//Jacobi et GS divergent

//question 3 
function X=Jacobi(A,b,X0,tol)
    M=diag(diag(A));
    N=M-A;
    incr=1;
    X=X0; Xaux=X0;
    while incr>tol
        X=M\(N*X+b);
        incr=norm(X-Xaux);
        Xaux=X;
    end
endfunction
A=[1 2 -2; 1 1 1; 2 2 1];
b=[1;1;1]
X=Jacobi(A,b,b,0.001)
disp("solution donnée par la méthode de Jacobi à 0.001 près")
disp(X)
disp("solution du systeme")
disp(A\b)