clear
lines(0)

//calcul le conditionnement de la matrice
//passe une matrice en un entier afin de pouvoir faire des comparaison
//Utilise la norme matricielle 2(norm())
function k=conditionnementMatrice(A)
    k=norm(A)*norm(inv(A))
endfunction

//Partie A

//Construit une matrice diagonal Delta Q°A.1
function D=matriceDiag(k,n)
    D=zeros(n,n)
    D(1,1)=1
    for i=2:n
        D(i,i)=k^(-1/(n-i+1))
    end
endfunction

//Construction de la matrice tridiagonal H Q°A.2.b
function H=matriceTriDiag(n)
    H=zeros(n,n)
    for i=1:n
        H(i,i)=2
    end
    for i=1:n-1
        H(i+1,i)=-1
        H(i,i+1)=-1
    end
endfunction

//Construction de la matrice A=U*Delta*V^T Q°A.2.d
function A=matriceAUDV(k, n)
    H=matriceTriDiag(n)
    [U,S,V]=svd(H)//U et V sont orthogonale
    D=matriceDiag(k, n)
    A=U*D*V'
endfunction

//Partie B

//Construit le second menbre b avec la solution exact epsilone[1,1...,1]^T Q°B.1.a (De la corection)
function b=secondMembre(A)
    n=size(A, 1) //permet de savoir les dimention d'epsilone
    epsilone=solutionEpsilone(n)
    b=A*epsilone
endfunction

//Construit epsilone[1,1...,1]^T Q°B.1.a (De la corection)
function epsilone=solutionEpsilone(n)
    epsilone=ones(n,1)
endfunction

//Resoud le SEL avec la bibliotheque Scilab (linsolve) Q°B.1.b
function x=resoudScilab(k,n)
    A=matriceAUDV(k,n)
    b=secondMembre(A)
    x=linsolve(A,-b)
endfunction

//Donne l'erreur de la solution Q°B.1.c
function e=calcErreur(k,n)
    x=resoudScilab(k,n)
    epsilone=solutionEpsilone(n)
    e=norm(epsilone-x)
endfunction

//Decompose une matrice en 3 Q°B.1.d
function [L,D,U]=decomposeLDU(A)
    [n,m]=size(A)
    L=zeros(n,m)
    D=zeros(n,m)
    U=zeros(n,m)
    for i=1:n
        for j=1:m
            if i>j then
                L(i,j)=A(i,j)
            elseif i==j then
                D(i,j)=A(i,j)
            else
                U(i,j)=A(i,j)
            end
        end
    end
endfunction

//Calcul l'erreur residuelle de la methode de jacobi Q°B.1.e
function r=erreurResidu(x,A,b)
    r=norm(A*x-b)/norm(b)
endfunction

//Resoud le SEL avec la methode de Jacobi Q°B.1.f
//i le nbrMax d'iteration et p la précision a atteindre définitive 
function x=resoudJacobi(k,n,i,p)
    A=matriceAUDV(k,n)
    b=secondMembre(A)
    x=zeros(n,1)//initialisation de la solution
    x=resoudJacobiContinu(x,A,b,i,p)
endfunction

//Pousuit les calculs si l'utilisateur le demande Q°B.1.h
//x est le x(k) calculer avant, i le nbrMax d'iteration, p la précision
function x=resoudJacobiContinu(x,A,b,i,p)
    j=i//sert si on doit répéter le calcul
    test=x//vérifie qu'il n'y a pas de divergence
    [L,D,U]=decomposeLDU(A)
    while j>0&erreurResidu(x,A,b)>p&erreurResidu(x,A,b)<=erreurResidu(test,A,b)
        //S'arrete si on dépasse i ou la précision demandée est atteinte
        test=x
        x=-inv(D)*(L+U)*x+inv(D)*b
        j=j-1
    end
    
    //Après calcul
    if erreurResidu(x,A,b)>erreurResidu(test,A,b) then//detecte la divergence
        k=conditionnementMatrice(A)
        n=size(A,1)
        disp('Il y a divergence pour k='+string(k)+' et n='+string(n))
        x=0
    else
        //x=postCalculJacobi(x,A,b,i,p,j)
    end
endfunction

//Affiche un compte rendu et propose si i=0 de continuer les calculs Q°B.1.h
function x=postCalculJacobi(x,A,b,i,p,j)
    eJ=erreurResidu(x,A,b)
    eD=erreurResidu(linsolve(A,-b),A,b)
    disp("Erreur résiduelle de la méthode jacobi est : "+string(eJ))
    disp('Erreur résiduelle de la méthode directe est : '+string(eD))
    disp("La différence est de : "+string(abs(eD-eJ)))
    disp("")
    if j<=0 then
        o=1//sert pour le input
        toc()//arrete le chrono pendant la saisi de l'utilisateur
        rep=input("La précision n''a pas été atteinte. Voulez vous continuer à calculer (o/n)?")
        tic()//reprend le chrono
        if rep==1 then
            x=resoudJacobiContinu(x,A,b,i,p)
        end
    end
endfunction

//Fonction principale qui genere les valeur k et n demande pour le SEL Q°B.1.h
function principale()
    k={10,10^3,10^8,10^16}
    n={10,50,100,150}
    for i=1:4
        for j=1:4
            resoudJacobi(k(i),n(j),300,10^-8)
        end
    end
endfunction

//Partie C

//Calcul la bonne valeur de m selon l'ennoncé Q°C
//H exposant
//On sait que ||e(0)||=1
//Si la valeur est négative ça veux dire qu'il y a divergence
function m=calculNbrIteration(k,n,H)
    Rj=calculRj(k,n,1)
    rayonSpec=calculRayonSpec(Rj)
    m=H/(log10(rayonSpec))
    m=ceil(m)
endfunction

//Calcul le Rj Q°C
//w (omega) sert à la Relaxation (par défault il faut 1)
function Rj=calculRj(k,n,w)
    A=matriceAUDV(k,n)
    [L,D,U]=decomposeLDU(A)
    Rj=(1-w)*eye(n)-w*inv(D)*(L+U)
endfunction

//Calcul le rayon spectral de A Q°C
function rayonSpec=calculRayonSpec(A)
    //lambda=spec(A)
    lambda=diag(bdiag(A))
    rayonSpec=max(abs(lambda(1)))
endfunction

//Calcul la bonne valeur de m selon l'ennoncé Q°D
function m=calculNbrIterationRelaxe(k,n,H)
    Rj=calculRjOptial(k,n)
    rayonSpec=calculRayonSpec(Rj)
    m=H/(log10(rayonSpec))
    m=ceil(m)
endfunction

//Calcul le Rj optimal en fonction de w Q°D
function RjOpt=calculRjOptial(k,n)
    w={0.5:0.1:1.8}
    RjOpt=calculRayonSpec(calculRj(k,n,w(1)))
    for i=2:size(w,2)
        if RjOpt>calculRayonSpec(calculRj(k,n,w(i))) then
            RjOpt=calculRayonSpec(calculRj(k,n,w(i)))
        end
    end
endfunction

//Fonction principale qui genere les valeur k et n demande pour le SEL Q°C et Q°D
function principale2()
    H=-8
    k={10,10^3,10^8,10^16}
    n={10,50,100,150}
    for i=1:4
        for j=1:4
            c='Pour k='+string(k(j))+'et n='+string(n(i))+' le nbr d''iteration '
            disp(c+'est de '+string(calculNbrIteration(k(j),n(i),H)))
            disp(c+'avec relaxation est de '+string(calculNbrIterationRelaxe(k(j),n(i),H)))
        end
    end
endfunction

//Donne un petit graphique pour le rapport (30s)
function graphique()
    k={10:10:500};
    n={10:5:150};
    z=zeros(size(k,2),size(n,2));
    for i=1:size(k,2)
        for j=1:size(n,2)
            z(i,j)=calculNbrIteration(k(i),n(j),-8);
        end
    end
    plot3d(k,n,z)
endfunction

function rep=divergence()
    i=300;
    p=10^-8;
    k={10,100};
    n={10:5:150};
    for i=1:size(k,2)
        for j=1:size(n,2)
            a=0==resoudJacobi(k(j),n(i),i,p)
            c='Pour k='+string(k(j))+'et n='+string(n(i))+' il y a '
            if a(1)then
                disp(c+' divergence')
            else
                disp(c+' convergence')
            end
        end
    end
    
endfunction




//Execution du programme principale
tic();

//principale()

//principale2()

//graphique()

k=10;
n=10;
i=300;
p=10^-8;

calculRjOptial(10,10)

divergence();




Information="Temps d''execution de l''algorithme : "+string(toc())+'s' //C'est plus joli dit ainsi

//La commande pour executer sans avoir tout le code affiché à chaque fois est : 
//exec('H:\Users\Administrator\Documents\Travail\ING-1\Analyse Numérique\TP Noté\Code\Code.sce', 0)
