n = 5;
A = rand(n,n); 

// Calcul des normes et du max_{i,j} abs(a_{ij}) 
m = max(A);
A1 = norm(A,1);
A2 = norm(A,2); // ou norm(A) 
Ainf = norm(A,'inf');
Afro = norm(A,'fro');
 
// On constate que m est inferieur a chacune de ces normes 
// Cette propriete se demontre facilement

////////////////////////////////////////////////////////////

// Fonction permettant le calcul du rayon spectral d'une matrice A. 

function [rs] = RayonSpectral(A)
 rs = norm(spec(A),'inf')
endfunction

function [resultat] = PuissanceNormep(A,p,k)
  resultat = norm(A^k,p)^(1./k)
endfunction

r_spec_A = RayonSpectral(A);

p = 1;
PuissanceNormep(A,p,1);
PuissanceNormep(A,p,2);
PuissanceNormep(A,p,10);
PuissanceNormep(A,p,20);
PuissanceNormep(A,p,50);
PuissanceNormep(A,p,70);
PuissanceNormep(A,p,100);
PuissanceNormep(A,p,200);
PuissanceNormep(A,p,500);

// On constate que PuissanceNormep(A,p,k) converge vers
// RayonSpectral(A) lorsque k tend vers l'infini. 

// Fonction permettant le calcul de l'inverse.
function B =Puissanceinverse(A)
  Sumpui = eye(A);
  C0 = eye(A) - A;
  C=C0;
  if (RayonSpectral(C) >= 1) then
    error( "le rayon spectral de (A-Id) > 1")
    else 
      while (norm(C) < 1000)&( norm(C) > 1.E-10) 
        Sumpui = Sumpui + C;
        C = C * C0;
      end;
      B=Sumpui 
  end
endfunction

// Matrices de Hilbert

function [H] = MatriceHilbert(n)
 H = zeros(n,n);
 for i=1:n
  for j=1:n
   H(i,j) = 1./(i+j-1);
  end
 end
endfunction

for k=1:10
 c(k) = cond(MatriceHilbert(k));
end

x = [1:10]';
xbasc(0);
plot2d([x,x],[log(c), 3.5*x-4.3]);

// On constate que le conditionnement en norme 2 de la 
// matrice de Hilbert d'ordre n est de l'ordre de exp(7n/2). 
// Cela se demontre rigoureusement... 

// Autre exemple de matrice mal conditionnee 

A1 = [10 7 8 7;7 5 6 5;8 6 10 9;7 5 9 10];
B1 = [25 -41 10 -6;-41 68 -17 10;10 -17 5 -3;-6 10 -3 2];
b1 = [32 23 33 31]';

x1 = linsolve(A1,-b1);

db = 0.01*rand(4,1);
//db = [0.1 -0.1 0.1 -0.1]';
b2 = b1 + db;

x2 = linsolve(A1,-b2);

[b1, db]
norm(db,2)/norm(b1,2)
[x1, x2-x1]
norm(x2-x1,2)/norm(x1,2)

dA = 0.01*rand(4,4);
A2 = A1 + dA;
x2 = linsolve(A2,-b1);

norm(dA,2)/norm(A1,2)
norm(x2-x1,2)/norm(x2,2)

// Verification des inegalites demontrees en cours

//norm(x2-x1,2)/norm(x1,2)
//cond(A1)*norm(db,2)/norm(b1,2)

//norm(x2-x1,2)/norm(x2,2)
//cond(A1)*norm(dA,2)/norm(A1,2)
 




