clear

stacksize('max')
format(25)

// A Construction des matrices 

function [Delta]=MatDiag(K,N)

Delta=zeros(N,N);
Delta(1,1)=1;
Delta(N,N)=1/K;

for i=2:N
    Delta(i,i)=K^(-1/(N-i+1))  // construit la matrice diagonale 
end

endfunction

function [K]=Veri(Delta)

K=norm(Delta,2)*norm(inv(Delta),2) // multiplie la norme 2 de A et la norme 2 de l'inverse de A
endfunction

function [H]=TriDiag(N)

H=diag(ones(1,N)*2)+diag(ones(1,N-1)*-1,-1)+diag(ones(1,N-1)*-1,1) // construit la matrice tri-diagonal 

endfunction

function [A]=Dec(H) // décomposition en valeur singulière 

[U,S,V]=svd(H)
A=U*Delta*V'

endfunction

function [Kb]=Verif(A)
    Kb=norm(A,2)*norm(inv(A),2) 
endfunction


K=10;
N=5;
delta= Mat(K,N)
H=TriM(N)
delta
H
[U,S,V]=svd(H)
A=U*delta*V'
A
k=Veri(A)
n=Veri(delta)

// B Calculs 

function [b]=ScdM(A,N)
    Q=ones(1,N)
    b=A*Q'
endfunction

function [x]=Solv(A,b)
    [x]=linsolve(A,-b)
endfunction

function [Er]=Erreur(N,x)
    Er=norm(ones(N,1)-(x))
endfunction
    
function [Jac]=Jacobi(k,A,N)
    
    D=diag(diag(A))
    L=tril(A)-D
    U=triu(A)-D
    b=A*ones(N,1)
     
     if k==0 then Jac=zeros(N,1)
        else Jac=-inv(D)*(L+U)*Jacobi(k-1,A,N)+inv(D)*b
    end
endfunction

function [i]=kJac(A,N,k)
     b=A*ones(N,1)
     i=0
     kJ=norm(A*Jacobi(i,A,N)-b,2)/norm(b,2)
     while kJ>10^-8 & k>i
         i=i+1 
         kJ=norm(A*Jacobi(i,A,N)-b,2)/norm(b,2)
     end
endfunction

function kJacobi=kJ(i,A,N)
     b=A*ones(N,1)
     kJacobi=norm(A*Jacobi(i,A,N)-b,2)/norm(b,2)    
endfunction

function [ErJ]=ErreurJ(i,A,N)
    ErJ=norm((ones(N,1))-(Jacobi(i,A,N)))
endfunction

function DErr=DErreur(ErJ,Er)
    DErr=ErJ/Er
endfunction

function [ibis]=kJacb(A,N,k)
     b=A*ones(N,1)
     ibis=0
     kJ=norm(A*Jacobi(ibis,A,N)-b,2)/norm(b,2)
     while kJ>10^-15 & k>ibis
         ibis=ibis+1 
         kJ=norm(A*Jacobi(ibis,A,N)-b,2)/norm(b,2)
     end
endfunction

b=ScdM(A,N)


// C. Erreur et nombre d'itérations

function [Lamb]=Spect(A)
    D=diag(diag(A))
    L=tril(A)-D
    U=triu(A)-D
    Rj=-inv(D)*(L+U)
    Lamb=abs(spec(Rj))
endfunction

function [LambM]=RSpec(Lamb)
    Lamb=gsort(Lamb)
    LambM=Lamb(1)    
endfunction

function [m]=mLamb(LambM)
    m=round(abs(8/log(LambM)))
endfunction

// D. Influence de la relaxation

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

function [Lambw]=Spectw(A,N,w)
    D=diag(diag(A))
    L=tril(A)-D
    U=triu(A)-D
    Rjw=(1-w)*diag(ones(1,N))-w*inv(D)*(L+U)
    Lambw=abs(spec(Rjw))
    Lambw=gsort(Lambw)
endfunction

function [LambM1]=RSpecw1(Lambw)
    LambM1=Lambw(N)    
endfunction

function [LambMN]=RSpecwN(Lambw)
    LambMN=Lambw(1)    
endfunction

function [w]=Val(LambM1,LambMN)
    w=2/(2-LambM1-LambMN)
endfunction