//Partie A

format(25);

//On créait une fonction afin de construire une matrice diagonale
function [Delta]=Diagonal(k,n)
    
    //on remplit la première valeure de la matrice
    D(1)=1;
    u=k
    v=n
        //on remplit la matrice
       for i=2:v
           D(i)=k^(-1/(n-1))
           n=n-1
       if i==u then
        D(i)=k^(-1)
           end
        end
    Delta=diag([D])
endfunction

//Phase de Test avec ici en exemple k = 2 et n = 6
k=10;
n=10;
D=Diagonal(k,n)

//On souhaite construire une matrice tridiagonale T avec m en paramétre

function [T]=tridiagonale(m)
    c=1;   
    //Rremplissage des différentes diagonales
    for i=1 :m do
      if (i==1) then
        T(i,c)=2;
        T(i,c+1)=-1;
      elseif (i==m) then
        T(i,c)=2;
        T(i,c-1)=-1;
      else
        T(i,c-1)=-1;
        T(i,c)=2;
        T(i,c+1)=-1;
       end
        c=c+1;
    end
endfunction

//Phase de Test en choisissant un m (exemple ici : m=10)
m=10;
t=tridiagonale(m);

// 
function [solution,erreur]=resolv(J,p)
    //C la matrice que l'on utilise dans l'équation
    //n la dimension de la matrice
    c=J*ones(p,1);
    solution=linsolve(J,-c);
    erreur=norm(solution-ones(p,1));
endfunction

// PARTIE B : Calculs à effectuer avec Scilab


function [erreur]=calculErreur(x,n,a)
    //on a x la solution précédente, n la dimension et a l'ordre
        erreur=norm((x-ones(n,1)),2);
endfunction

//on créait la fonction Jacobi pour utiliser la méthode itérative de Jacobi

function [fctJacobi]=Jacobi(k,A,n)
    //Avec k la valeurs
    //A la matrice de départ et n sa dimension
    
    D=diag(diag(A))
    L=tril(A)-D
    U=triu(A)-D
    b=A*ones(n,1)
     
     if k==0 then fctJacobi=zeros(n,1)
        else fctJacobi=-inv(D)*(L+U)*fctJacobi(k-1,A,n)+inv(D)*b
    end
endfunction

//L'itération de Jacobi nous permet de calculer l'erreur raltive avec une certaine précision
function [i]=iterationJacobi(A,n,k)
    //k le nombre d'itération voulu
     b=A*ones(n,1)
     i=0
     errResiduelle=norm(A*Jacobi(i,A,n)-b,2)/norm(b,2)
     while kJ>10^-8 & k>i
         i=i+1 
         errResiduelle=norm(A*Jacobi(i,A,n)-b,2)/norm(b,2)
     end
endfunction


// A présent on calcul l'erreur relative de la solution
function [kJacobi]=errResiduelle(i,A,n)
    //i la valeur que l'on souhaite
     b=A*ones(n,1)
     kJacobi=norm(A*Jacobi(i,A,n)-b,2)/norm(b,2)    
endfunction

//Calcul de l'erreur suite à la méthode Jacobi
function [ErJ]=ErreurJacobi(i,A,N)
    ErJ=norm((ones(N,1))-(Jacobi(i,A,N)))
endfunction


//
function [erreurRelative]=erreur(ErJacobi,Erreur)
    erreurRelative=/Erreur
endfunction


//par cette fonction on continue les calculs pour améliorer l'erreur résiduelle qui va se rapprocher du seuil
function [ibis]=iterationJacobiBis(A,n,k)

     b=A*ones(n,1)
     ibis=0
     errResiduellebis=norm(A*Jacobi(ibis,A,n)-b,2)/norm(b,2)
     while kJ>10^-8 & k>i
         ibis=ibis+1 
         errResiduellebis=norm(A*Jacobi(ibis,A,n)-b,2)/norm(b,2)
     end
endfunction

// PARTIE C : Erreur et nombre d'itérations

//Obtention des Valeurs propres (en val absolue) de Rj matrice d'itération de Jacobi
function [Lambda]=Spectre(A)
    D=diag(diag(A));
    L=tril(A)-D;
    U=triu(A)-D;
    Rj=-inv(D)*(L+U);
    Lambda=abs(spec(Rj));
endfunction


//Trie des valeurs propres dans l'ordre décroissant
//lambda est le vecteur créait par ce trie
//LambdaMax est le max de ces valeurs propres
function [LambdaMax]=RayonSpectral(Lambda)

    Lambda=gsort(Lambda);
    LambdaMax=Lambda(1);
endfunction


//on cherche le nombre d'itération m nécessaire pour faire diminuer l'erreur
function [m]=mLambda(LambMax)
    m=round(abs(8/log(LambMax)))
endfunction

// PARTIE D : Influence de la relaxation


//Cette fonction nous permet d'obtenir les valeurs des valeurs propres de Rj en valeur absolu ou Rj est la matrice d'itération de Jacobi
function [Jacw]=Jacobiw(k,A,N,w)
    
    D=diag(diag(A));
    L=tril(A)-D;
    U=triu(A)-D;
    b=A*ones(N,1);
     
     if k==0 then Jacw=zeros(N,1);
        else Jacw=[(1-w)*diag(ones(1,N))-w*inv(D)*(L+U)]*Jacobiw(k-1,A,N,w)+w*inv(D)*b;
    end
endfunction

    //cette fonction nous donne un vecteur trié dans l'ordre decroissant des valeurs propres de Rj avec le coefficient de relaxation utilisé lors du calcul.
function [Lambw]=Spectw(A,n,w)
    D=diag(diag(A))
    L=tril(A)-D
    U=triu(A)-D
    Rjw=(1-w)*eye(1,n))-w*inv(D)*(L+U)
    Lambw=abs(spec(Rjw))
    Lambw=gsort(Lambw)
endfunction


    //cette fonction renvoie la plus petite valeur propre du vecteur Lambdaw qui a été trié précedemment
function [LambdaMin1]=minSpectrew1(Lambdaw)
    LambdaMin1=Lambdaw(N)    
endfunction



function [LambdaMaxN]=maxSpectrewN(Lambdaw)
        //cette fonction renvoie la plus grande valeur propre du vecteur Lambdaw qui a été trié précedemment
    //lambdaw le vecteur contenant les valeurs propres triées dans l'ordre decroissant
    LambdaMaxN=Lambdaw(1)    
endfunction


    //cette fonction va nous permettre de comparer notre valeur de relaxation choisit au départ avec les valeurs calculés via les lambdas.
    //w la relaxation
function [w]=Validite(LambdaMin,LambdaMax)
    w=2/(2-LambdaMin-LambdaMax)
endfunction