/////////////////////////////////Partie A/////////////////////////////////// Attention les gars changer le code!! Il manque la partie D


//Fonction permettant de construire la matrice diagonale

function [Delta]=Diagonal(k,n)
L(1)=1;
u=k
v=n
    for i=2:v
       L(i)=k^(-1/(n-1))
       n=n-1
       if i==u then
           L(i)=k^(-1)
       end
    end
    Delta=diag([L])
endfunction

//Fonction permettant d'obtenir la matrice tridiagonale avec comme paramètre n qui correspond à la taille de la matrice
function [H]=tridiagonale(n)
    c=1;  
    for i=1 :n do
            if (i==1) then
                H(i,c)=2;
                H(i,c+1)=-1;
            elseif (i==n) then
                H(i,c)=2;
                H(i,c-1)=-1;
            else
                H(i,c-1)=-1;
                H(i,c)=2;
                H(i,c+1)=-1;
            end
            c=c+1;
    end
endfunction


// On test comme suit dans la console les valeurs suivantes :
//Pour la fonction diagonale on a utilisé k=10 et n=10 puis on a stocké H fonction tridiagonale où on a vérifié les calculs demandés

printf("\n------------------------Partie A: Matrice à conditionnement donné--------------\n");
k=10;//valeur du conditionnement, cette valeurs peut être modifier, on va la faire passer 10000 dans la partie D
n=10;//taille de la matrice, cette valeurs peut être modifier, on va la faire passer 150 dans la partie D
printf("\nOn choisit n= %d et k= %d ",n,k);
printf("\nConditionnement de la matrice A:%d \n",k);


Delta=Diagonal(k,n);
H=tridiagonale(n);
[U,S,V]=svd(H);
A=U*Delta*V';

//////////////////////////////////Partie B/////////////////////////////////

function [b]=segondmembre(n)
    for (i=1:n)
        SE(i,1)=1
    end
    b=A*SE
endfunction

b=segondmembre(n)
//-----------------méthode directe + erreurs----------------

[x]=linsolve(A,-b);//résolution du système avec la méthode directe
erreurD=norm(ones(n,1)-x,2);// erreur de la méthode directe
erreurResD=norm(b-A*x,2); // erreur résiduelle de la méthode directe

printf("\n---------------------Partie B: Erreur résiduelle de la méthode directe et itérative----------------------\n");
printf("\nErreur résiduelle direct: %.20f",erreurResD);

//-----------------méthode jacobienne +erreurs------------------

D=diag(diag(A))//la diagonale de la matrice A
L=tril(A)-D//la partie de A strictement inférieur à la diagonale
U=triu(A)-D//la partie de A inférieure à la diagonale
i=1000//le nombre d'itération utilisé dans la fonction des erreurs. On ne pouvais pas prendre 200 000 comme dit dans l'énoncé.

//La méthode itérative de Jacobi
function [X]=Jacobi(i)
    X0=zeros(n,1);
    X1=(D^-1)*b;
    res=(-D^-1)*(L+U)*X0+X1;
   
    if (i==0) then
        X=X1;
    else
        for i=1:(i-1)
            res=(-D^-1)*(L+U)*res+(X1);
        end
        X=res
    end
endfunction

X=Jacobi(k);

function [erreur]=rk(i)//erreur résiduelle de la méthode Jacobi
    k=0;
   erreur=(norm(A*Jacobi(k)-b,2)/norm(b,2));
    while erreur>10^-17 & i>k //on remarque que les valeurs à partir de 10^-16 restent inchangé
        k=k+1;
        erreur=(norm(A*Jacobi(k)-b,2)/norm(b,2));
    end
endfunction



erreurResJ=rk(i);

printf("\nErreur résiduelle avec %d itérations : %.20f",i,erreurResJ);
printf("\nComparer les deux erreurs résiduelles en faisant le rapport des deux erreurResJ/erreurResD: %f\n",erreurResJ/erreurResD);
printf("on a donc %f<1, donc erreurResJ<erreurResD\n",erreurResJ/erreurResD);


function [erreurJ]=ek(k)// erreur de la méthode jacobienne
   erreurJ=norm(ones(n,1)-Jacobi(k),2);
endfunction

erreurJ=ek(k)


//////////////////////////////////////////////////Partie C///////////////////////////////////////////////////

function [h]=RechercheH()
    global('lambda');
    global('rayonspectral');
     Rj= (-D^-1)*(L+U);  //Matrice d'itération de Jacobi
    [lambda]= diag(bdiag(Rj));//calcul du rayon spectral de la matrice Rj
    rayonspectral= max(abs(lambda));
    htemp=1;
    m=norm(htemp/log(rayonspectral),2);
    while ek(m)>10^-8
        global('m');
        m=norm(htemp/log(rayonspectral),2);
        htemp=htemp+1;
    end
    h=htemp;
endfunction



printf("\n------------------Partie C: Erreur et nombre ditération---------------------\n");
global('m');
global('rayonspectral');
Hj=RechercheH();

printf("\non obtient H= %d et m= %d itérations pour que Ek soit < 10^-8\n",Hj,m);