//////////////////////////////
// Auteurs : Vanhove Romain
//           Bertiaux Mathias
// Ing1 groupe B14
//

clear();
format(25);
stacksize('max')
 
 
 
 // Fonction qui crée la matrice diagonale
 
  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));
      end
  endfunction
  
  
  

  

// Fonction qui calcule le conditionnement de la matrice A

  function [K]=Conditionnement(A)
      K=norm(A,2)*norm(inv(A),2);
  endfunction
  
  
  

  

// Fonction qui construit la matrice tridiagonale

  function [H]=TriDiag(N)
      H=diag(ones(1,N)*2)+diag(ones(1,N-1)*-1,-1)+diag(ones(1,N-1)*-1,1);
  endfunction
  
  
  

  

// Fonction qui fait une décomposition en valeurs singulières

  function [A]=Decomposition(H)
      [U,S,V]=svd(H);
      A=U'*delta*V; // delta définit plus haut
   endfunction
  
  
  

  

//Fonction qui calcule b=A*x avec x = [1,1,1, ... , 1, 1]'

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

  
   
   // Fonction qui résoud le système A*x = b
   
  function [x]=Sol(A,b)
       x=linsolve(A,-b);
  endfunction
   
  
  

  
 
 // Fonction qui calcule l'erreur, première méthode
 
  function [E]=Erreur(N,A,b) 
      tmp=Sol(A,b);
      E=norm(ones(N,1)-(tmp));
  endfunction
   
  
  

  

 // Fonction qui calcule l'erreur, deuxième méthode
 
  function [E]=Erreur2(A,b) 
      [x, fail, E, iter, res] = pcg(A, b, 10^(-12), 10);
  endfunction
   
  
  

  

  // Fonction qui calcule l'erreur, troisième et quatrième méthode (il suffit de changer la matrice T)
 
  function [E]=Erreur3(A,b, T) 
      [x, fail, E, iter, res] = pcg(A, b, 10^(-12), 10, T);
  endfunction
   
  
  

  

 // Fonction main
    
    function [E4, E3, E2, E1]=main(K, N)
        delta=MatDiag(K,N);
        H=TriDiag(N);
        A=Decomposition(H);
        b=Calculb(A,N);
        //Premiere méthode
        E1=Erreur(N,A,b)
        //Deuxième méthode    
        E2=Erreur2(A,b)
        //Troisième méthode      
        P1=diag(diag(A));
        E3=Erreur3(A,b,P1)
        //Quatrième méthode      
        D=diag(diag(A));
        L=tril(A,-1);
        P2=(D-L)*inv(D)*(D-L');
        E4=Erreur3(A,b,P2)
    endfunction
      
     	disp("============================================================================================")

	      disp("Les erreurs relatives à chaque méthode de résolution sont listées de la manière suivante :")
	      disp("E1, erreur de la méthode de résolution directe avec linsolve,")
	      disp("E2, erreur de la méthode du gradient conjugué sans pré-conditionnement,")
	      disp("E3, erreur de la méthode du gradient conjugué avec pré-conditionnement, avec P=diag(A),")
	      disp("E4, erreur de la méthode du gradient conjugué avec pré-conditionnement, avec P=(D-L)*inv(D)*(D-L^t).")

	disp("============================================================================================")
	disp("")
      
  
	tmp1=input("Dimension de la matrice ? ");
	tmp2=input("Valeur de conditionnement ? ");
	disp("")


   	[E4, E3, E2, E1]=main(tmp2,tmp1)


 
