Chapitre 14  Algorithmes d’algébre

14.1  Méthode pour résoudre des systèmes linéaires

Dans Xcas, il existe déjà les fonctions qui correspondent aux algorithmes qui suivent, ce sont : ref, rref, ker, pivot

14.1.1  Le pivot de Gauss quand A est de rang maximum

Étant donné un système d’équations noté AX=b, on cherche à le remplacer par un système équivalent et triangulaire inférieur.
À un système d’équations AX=b, on associe la matrice M formée par A que l’on borde avec b comme dernière colonne.
La méthode de Gauss (resp Gauss-Jordan) consiste à multiplier A et b (donc M) par des matrices inversibles, afin de rendre A triangulaire inférieure (resp diagonale). Cette transformation, qui se fera au coup par coup en traitant toutes les colonnes de A (donc toutes les colonnes de M, sauf la dernière), consiste par des combinaisons de lignes de M à mettre des zéros sous (resp de part et d’autre) la diagonale de A.
La fonction gauss_redi ci-dessous, transforme M par la méthode de Gauss, la variable pivo (car pivot est une fonction de Xcas) sert à mettre le pivot choisi.
Pour j=0..p−2 (p−2 car on ne traite pas la dernière colonne de M), dans chaque colonne j, on cherche ce qui va faire office de pivot à partir de la diagonale : dans le programme ci-dessous on choisit le premier élément non nul, puis par un échange de lignes, on met le pivot sur la diagonale (pivo=M[j,j]). Il ne reste plus qu’à former, pour chaque ligne k (k>j) et pour a=M[k,j], la combinaison :
pivo*ligne_k−a*ligne_j (et ainsi M[k,j] devient nul).
On écrit pour réaliser cette combinaison :
a:=M[k,j];
for (l:=j;l<nc;l++){
M[k,l]:=M[k,l]*pivo-M[j,l]*a;}
On remarquera qu’il suffit que l parte de j car :
pour tout l<j, on a déjà obtenu, par le traitement des colonnes l=0..j−1, M[k,l]=0.
Le programme ci-dessous ne sera utile que si on trouve un pivot pour chaque colonne : c’est à dire si la matrice A est de rang maximum. En effet, dans ce programme, si on ne trouve pas de pivot (i. e. si tous les éléments d’une colonne sont nuls sur et sous la diagonale), on continue comme si de rien était...

gauss_redi(M):={
local pivo,j,k,nl,nc,temp,l,n,a;
nl:=nrows(M);
nc:=ncols(M);
n:=min(nl,nc-1);
//on met des 0 sous la diagonale du carre n*n
for (j:=0;j<n;j++) {
  //choix du pivot mis ds pivo
    k:=j;
    while (M[k,j]==0 and k<nl-1) {k:=k+1;}
    //on ne fait la suite que si on a pivo!=0
    if (M[k,j]!=0) {
      pivo:=M[k,j];
      //echange de la ligne j et de la ligne k
      for (l:=j;l<nc;l++){
       temp:=M[j,l];
       M[j,l] := M[k,l];
        M[k,l]:=temp;     
      }
     //fin du choix du pivot qui est M[j,j]
      for (k:=j+1;k<nl;k++) {
        a:=M[k,j];
        for (l:=j;l<nc;l++){ 
          M[k,l]:=M[k,l]*pivo-M[j,l]*a; 
         } 
       }
    }  
}
return M;
}

On tape :
M0:= [[1,2,3,6],[2,3,1,6],[3,2,1,6]]
gauss_redi(M0)
On obtient :
[[1,2,3,6],[0,-1,-5,-6],[0,0,-12,-12]]
On tape :
M1:= [[1,2,3,4],[0,0,1,2],[0,0,5,1]]
gauss_redi(M1)
On obtient :
[[1,2,3,4],[0,0,1,2],[0,0,5,1]]
On tape :
M2:= [[1,2,3,4],[0,0,1,2],[0,0,5,1],[0,0,3,2],[0,0,-1,1]]
gauss_redi(M2) On obtient :
[[1,2,3,4],[0,0,1,2],[0,0,5,1],[0,0,0,7],[0,0,0,6]]
c’est à dire :

14.1.2  Le pivot de Gauss pour A quelconque

On cherche un programme valable lorsque A est quelconque : on veut, dans ce cas, mettre des zéros sous la "fausse diagonale" de A (ce qu’on appelle "fausse diagonale" c’est la diagonale obtenue en ne tenant pas compte des colonnes pour lesquelles la recherche du pivot a été infructueuse : un peu comme si ces colonnes étaient rejetées à la fin de la matrice).
Donc, dans ce programme, si on ne trouve pas de pivot pour la colonne d’indice jc (i. e. si tous les éléments de la colonne jc sont nuls sur et sous la diagonale), on continue en cherchant, pour la colonne suivante (celle d’indice jc+1), un pivot à partir de l’élément situé à la ligne d’indice jc (et non comme précédemment à partir de jc+1), pour mettre sur la colonne jc+1, des zéros sur les lignes jc+1,...,nl−1. On est donc obligé, d’avoir 2 indices jl et jc pour repérer les indices de la "fausse diagonale".

gauss_red(M):={
local pivo,jc,jl,k,nl,nc,temp,l,a;
nl:=nrows(M);
nc:=ncols(M);
//on met des 0 sous la fausse diagonale d'indice jl,jc 
jc:=0;
jl:=0;
//on traite chaque colonne (indice jc)
while (jc<nc-1 and jl<nl-1) {
  //choix du pivot que l'on veut mettre en M[jl,jc]
    k:=jl;
    while (M[k,jc]==0 and k<nl-1) {k:=k+1;}
    //on ne fait la suite que si M[k,jc](=pivo)!=0
    if (M[k,jc]!=0) {
      pivo:=M[k,jc];
      //echange de la ligne jl et de la ligne k
      for (l:=jc;l<nc;l++){
       temp:=M[jl,l];
       M[jl,l] := M[k,l];
        M[k,l]:=temp;     
      }
     //fin du choix du pivot qui est M[jl,jc]
      //on met des 0 sous la fausse diagonale de 
      //la colonne jc
      for (k:=jl+1;k<nl;k++) {
        a:=M[k,jc];
        for (l:=jc;l<nc;l++){ 
          M[k,l]:=M[k,l]*pivo-M[jl,l]*a; 
         } 
       }
      //on a 1 pivot,l'indice-ligne de 
      //la fausse diag augmente de 1
      jl:=jl+1;
    }//fin du if (M[k,jc]!=0)
    //colonne suivante,l'indice-colonne de
    //la fausse diag augmente de 1
    jc:=jc+1;  
}//fin du while
return M;
}

On tape :
M0:= [[1,2,3,6],[2,3,1,6],[3,2,1,6]]
gauss_red(M0)
On obtient :
[[1,2,3,6],[0,-1,-5,-6],[0,0,-12,-12]]
On tape :
M1:= [[1,2,3,4],[0,0,1,2],[0,0,5,1]]
gauss_red(M1)
On obtient :
[[1,2,3,4],[0,0,1,2],[0,0,0,-9]]
On tape :
M2:= [[1,2,3,4],[0,0,1,2],[0,0,5,1],[0,0,3,2],[0,0,-1,1]]
gauss_red(M2)
On obtient :
[[1,2,3,4],[0,0,1,2],[0,0,0,-9],[0,0,0,-4],[0,0,0,3]]

14.1.3  La méthode de Gauss-Jordan

Étant donnée un système d’équations, on cherche à le remplacer par un système diagonale équivalent. À un système d’équations AX=b, on associe la matrice M formée de A, bordée par b comme dernière colonne.
La fonction gaussjordan_redi transforme M par la méthode de Gauss-Jordan. On cherche dans chaque colonne j,(j=0..nc−2) à partir de la diagonale, ce qui va faire office de pivot : dans le programme ci-dessous on choisit le premier élément non nul, puis par un échange de lignes, on met le pivot sur la diagonale (pivo=M[j,j]). Il ne reste plus qu’à former, pour chaque ligne k (k≠ j) et pour a=M[k,j], la combinaison :
pivo*ligne_k−a*ligne_j (et ainsi M[k,j] devient nul).
On écrit pour réaliser cette combinaison :
lorsque k<j
a:=M[k,j];
for (l:=0;l<j;l++){
M[k,l]:=M[k,l]*pivo-M[j,l]*a;}
et lorsque k>j
a:=M[k,j];
for (l:=j;l<nc;l++){
M[k,l]:=M[k,l]*pivo-M[j,l]*a;}
On remarquera que l part soit de 0 soit de j car pour l<j, on a M[k,l]=0 seulement si k>j.
Si on ne trouve pas de pivot, on continue comme si de rien n’était : on obtiendra donc des zéros au dessus de la diagonale que si on a trouvé un pivot pour chaque colonne.

gaussjordan_redi(M):={
local pivo,j,k,nl,nc,temp,l,n,a;
nl:=nrows(M);
nc:=ncols(M);
n:=min(nl,nc-1);
//on met 0 sous et au dessus de la diag du carre n*n
for (j:=0;j<n;j++) {
  //choix du pivot
    k:=j;
    while (M[k,j]==0 and k<nl-1) {k:=k+1;}
    //on ne fait la suite que si on a pivo!=0
    if (M[k,j]!=0) {
      pivo:=M[k,j];
      //echange de la ligne j et de la ligne k
      for (l:=j;l<nc;l++){
       temp:=M[j,l];
       M[j,l] := M[k,l];
        M[k,l]:=temp;     
      }
     //fin du choix du pivot qui est M[j,j]
     // on met des zeros au dessus de la diag
     for (k:=0;k<j;k++) {
        a:=M[k,j];
        for (l:=0;l<nc;l++){ 
          M[k,l]:=M[k,l]*pivo-M[j,l]*a; 
         } 
      }
      // on met des zeros au dessous de la diag
      for (k:=j+1;k<nl;k++) {
        a:=M[k,j];
        for (l:=j;l<nc;l++){ 
          M[k,l]:=M[k,l]*pivo-M[j,l]*a; 
         } 
       }
    }  
}
return M;
}

De la même façon que pour la mèthode de Gauss, on va mettre des zéros sous la "fausse diagonale" et au dessus de cette "fausse diagonale" (on ne pourra pas bien sûr, mettre des zéros au dessus de cette "fausse diagonale", pour les colonnes sans pivot!!)

gaussjordan_red(M):={
local pivo,jc,jl,k,nl,nc,temp,l,a;
nl:=nrows(M);
nc:=ncols(M);
//on met des 0 sous la fausse diagonale 
jc:=0;
jl:=0;
//on doit traiter toutes les colonnes sauf la derniere 
//on doit traiter toutes les lignes
while (jc<nc-1 and jl<nl) {
  //choix du pivot que l'on veut mettre en M[jl,jc]
    k:=jl;
    while (M[k,jc]==0 and k<nl-1) {k:=k+1;}
    //on ne fait la suite que si on a pivo!=0
    if (M[k,jc]!=0) {
      pivo:=M[k,jc];
      //echange de la ligne jl et de la ligne k
      for (l:=jc;l<nc;l++){
       temp:=M[jl,l];
       M[jl,l] := M[k,l];
        M[k,l]:=temp;     
      }
     //fin du choix du pivot qui est M[jl,jc]
      //on met des 0 au dessus la fausse diagonale de 
      //la colonne jc
      for (k:=0;k<jl;k++) {
        a:=M[k,jc];  
        for (l:=0;l<nc;l++){ 
          M[k,l]:=M[k,l]*pivo-M[jl,l]*a; 
         } 
       }
       //on met 0 sous la fausse diag de la colonne jc
       for (k:=jl+1;k<nl;k++) {
        a:=M[k,jc];
        for (l:=jc;l<nc;l++){ 
          M[k,l]:=M[k,l]*pivo-M[jl,l]*a; 
         } 
       }
      //on a un pivot donc le numero de  
      //la ligne de la fausse diag augmente de 1
      jl:=jl+1;
    }
    //ds tous les cas, le numero de
    //la colonne de la fausse diag augmente de 1
    jc:=jc+1;  
}
return M;
}

On tape :
M0:= [[1,2,3,6],[2,3,1,6],[3,2,1,6]]
gaussjordan_red(M0)
On obtient :
[[12,0,0,12],[0,12,0,12],[0,0,-12,-12]]
On tape :
M1:= [[1,2,3,4],[0,0,1,2],[0,0,5,1]]
gaussjordan_red(M1)
On obtient :
[[1,2,0,-2],[0,0,1,2],[0,0,0,-9]]
On tape :
M2:= [[1,2,3,4],[0,0,1,2],[0,0,5,1],[0,0,3,2],[0,0,-1,1]]
gaussjordan_red(M2)
On obtient :
[[1,2,0,-2],[0,0,1,2],[0,0,0,-9],[0,0,0,-4],[0,0,0,3]]

14.1.4  La méthode de Gauss et de Gauss-Jordan avec recherche du pivot

On peut bien sûr modifier la recherche du pivot.
Pour les méthodes numériques il est recommandé de normaliser les équations (on divise chaque ligne k par max_j|a_k,j|) et de choisir le pivot qui a la plus grande valeur absolue.
En calcul formel, on prend l’expression exacte la plus simple possible. Ici on normalise les équations à chaque étape et on choisit le pivot qui a la plus grande valeur absolue : c’est ce que font, avec ce choix du pivot les programmes (cf le répertoire exemples), gauss_reducni (analogue à gauss_redi), gauss_reducn (analogue à gauss_red)), gaussjordan_reducni (analogue à gaussjordan_redi) et gaussjordan_reducn (analogue à gaussjordan_red) :

     
gaussjordan_reducn(M):={
local pivo,j,jc,jl,k,nl,nc,temp,l,a,piv,kpiv,maxi;
nl:=nrows(M);
nc:=ncols(M);
//on met des 0 sous la fausse diagonale 
jc:=0;
jl:=0;
while (jc<nc-1 and jl<nl) {
   //on normalise les lignes
   for (jj:=0;jj<nl;jj++) {
       maxi:=max(abs(seq(M[jj,kk], kk=0..nc-1)));
       for (kk:=0;kk<nc;kk++) {
           M[jj,kk]:=M[jj,kk]/maxi;
       }  
   }
   //choix du pivot que l'on veut mettre en M[jl,jc]
    kpiv:=jl;
    piv:=abs(M[kpiv,jc]);
    for (k:=jl+1;k<nl;k++){
    if (abs(M[k,jc])>piv) {piv:=abs(M[k,jc]);kpiv:=k;}
    }
   //on ne fait la suite que si on a piv!=0
    if (piv!=0) {
      pivo:=M[kpiv,jc];
      k:=kpiv;
      //echange de la ligne jl et de la ligne k
      for (l:=jc;l<nc;l++){
       temp:=M[jl,l];
       M[jl,l] := M[k,l];
        M[k,l]:=temp;     
      }
     //fin du choix du pivot qui est M[jl,jc]
     //on met des 0 au dessus la fausse diagonale 
     //de la colonne jc
      for (k:=0;k<jl;k++) {
        a:=M[k,jc];  
        for (l:=0;l<nc;l++){ 
          M[k,l]:=M[k,l]*pivo-M[jl,l]*a; 
         } 
       }
       //on met des 0 sous la fausse dia de la colonne jc
       for (k:=jl+1;k<nl;k++) {
        a:=M[k,jc];
        for (l:=jc;l<nc;l++){ 
          M[k,l]:=M[k,l]*pivo-M[jl,l]*a; 
         } 
       }
      //on a un pivot donc,le numero de
      //la ligne de la fausse diag augmente de 1
      jl:=jl+1;
    }
    //ds tous les cas,le numero de
    //la colonne de la fausse diag augmente de 1
    jc:=jc+1;  
}
return M;
}

On tape :
M0:= [[1,2,3,6],[2,3,1,6],[3,2,1,6]]
gaussjordan_reducn(M0)
On obtient :
[[5/6,0,0,5/6],[0,5/6,0,5/6],[0,0,1,1]]
On tape :
M1:= [[1,2,3,4],[0,0,1,2],[0,0,5,1]]
gaussjordan_reducn(M1)
On obtient :
[[1/4,1/2,0,17/20],[0,0,1,1/5],[0,0,0,9/10]]
On tape :
M2:=[[1,2,3,4],[0,0,1,2],[0,0,5,1],[0,0,3,2],[0,0,-1,1]]
gaussjordan_reducn(M2)
On obtient :
[[1/4,1/2,0,17/20],[0,0,1,1/5],
[0,0,0,9/10],[0,0,0,7/15],[0,0,0,6/5]]

On a donc :

14.1.5  Application : recherche du noyau grâce à Gauss-Jordan

Soit une application linéaire f de ℝⁿ dans ℝ^p, de matrice A dans la base canonique de ℝⁿ et de ℝ^p.
Trouver le noyau de f revient à résoudre f(X)=0 ou encore à résoudre AX=0 pour X∈ℝⁿ.
Pour cela, on va utiliser la méthode de Gauss-Jordan en rajoutant des lignes de 0 aux endroits où l’on n’a pas trouvé de pivot de façon à ce que les pivots suivants se trouve sur la diagonale (et non sur la "fausse diagonale").
Plus précisément :
- quand on a trouvé un pivot pour une colonne, à l’aide de la méthode habituelle de réduction de Gauss-Jordan, on met sur cette colonne :
1 sur la diagonale et
0 de part et d’autre du 1.
- quand on n’a pas trouvé de pivot pour la colonne j, on rajoute une ligne de 0 : la ligne j devient une ligne de 0 et les autres lignes sont décalées.
- on rajoute éventuellement à la fin, des lignes de 0, pour traiter toutes les colonnes de A et avoir une matrice carrée.
- on enlève éventuellement à la fin, des lignes de 0, pour avoir une matrice carrée.
Remarque : on ne change pas le nombre de colonnes.
Ainsi, si la colonne C_j a un 0 sur sa diagonale, si e_j est le j-ième vecteur de la base canonique de ℝⁿ, alors N_j=C_j−e_j est un vecteur du noyau. En effet, on a A(e_j)=C_j; si on a C_j=a1e₁+....+a_j−1e_j−1, pour k<j, et pour a_k!=0, on a A(e_k)=e_k c’est à dire A(C_j)=C_j=A(e_j), soit A(C_j−e_j)=A(N_j)=0.
Voici le programme correspondant en conservant la variable jl afin de faciliter la lisibilité du programme :

     
gaussjordan_noyau(M):={
local pivo,jc,jl,k,j,nl,nc,temp,l,a,noyau;
nl:=nrows(M);
nc:=ncols(M);
//on met des 0 sous la diagonale 
jc:=0;
jl:=0;
// on traite toutes les colonnes
while (jc<nc and jl<nl) {
  //choix du pivot que l'on veut mettre en M[jl,jc]
    k:=jl;
    while (M[k,jc]==0 and k<nl-1) {k:=k+1;}
    //on ne fait la suite que si on a pivo!=0
    if (M[k,jc]!=0) {
      pivo:=M[k,jc];
      //echange de la ligne jl et de la ligne k
      for (l:=jc;l<nc;l++){
       temp:=M[jl,l];
       M[jl,l] := M[k,l];
        M[k,l]:=temp;     
      }
     //fin du choix du pivot qui est M[jl,jc]
      //on met 1 sur la diagonale de la colonne jc
      for (l:=0;l<nc;l++) {
      M[jl,l]:=M[jl,l]/pivo;
      }
      //on met des 0 au dessus de la diagonale
      // de la colonne jc
      for (k:=0;k<jl;k++) {
        a:=M[k,jc];  
        for (l:=0;l<nc;l++){
          M[k,l]:=M[k,l]-M[jl,l]*a; 
         } 
       }
       //on met des 0 sous la diag de la colonne jc
       for (k:=jl+1;k<nl;k++) {
        a:=M[k,jc];
        for (l:=jc;l<nc;l++){ 
          M[k,l]:=M[k,l]-M[jl,l]*a; 
         } 
       }
    }
    else{
    //on ajoute une ligne de 0 si ce n'est pas le dernier 0
    if (jl<nc-1){
     for (j:=nl;j>jl;j--){
      M[j]:=M[j-1];
      }
    M[jl]:=makelist(0,1,nc);
    nl:=nl+1;
    }
   }
    //ds tous les cas,le numero de colonne et 
    //le numero de ligne augmente de 1
    jc:=jc+1;  jl:=jl+1;
    //il faut faire toutes les colonnes 
    if (jl==nl and jl<nc) { M[nl]:=makelist(0,1,nc);nl:=nl+1;}
}
noyau:=[];
//on enleve les lignes en trop pour avoir 
//une matrice carree de dim nc
//on retranche la matrice identite
M:=M[0..nc-1]-idn(nc);
for(j:=0;j<nc;j++){
if (M[j,j]==-1) {noyau:=append(noyau,M[0..nc-1,j]);}
}
return noyau;
}

Remarque
On peut écrire le même programme en supprimant la variable jl puisque jl=jc (on met jc à la place de jl et on supprime jl:=0 et jl:=jl+1).
On met ce programme dans un niveau éditeur de programmes (que l’on ouvre avec Alt+p), puis on le teste et on le valide avec OK.
On tape :
gaussjordan_noyau([[1,2,3],[1,3,6],[2,5,9]])
On obtient :
[[-3,3,-1]]
On tape :
gaussjordan_noyau([[1,2,3,4],[1,3,6,6],[2,5,9,10]])
On obtient :
[[-3,3,-1,0],[0,2,0,-1]]

14.2  Résolution d’un système linéaire

14.2.1  Résolution d’un système d’équations linéaires

L’algorithme

On associe à un système d’équations linéaires, une matrice A constituèe de la matrice du système augmentèe d’une colonne formée par l’opposé du second membre.
Par exemple au système : [x+y=4,x−y=2] d’inconnues [x,y] on associe la matrice :
A=[[1,1,−4],[1,−1,−2]].
Puis on réduit avec la méthode de Gauss-Jordan la matrice A pour obtenir une matrice B. Pour chaque ligne de B :
- si il n’y a que des zéro on regarde la ligne suivante,
- si il y a des zéro sauf en dernière position, il n’y a pas de solution,
- dans les autres cas on obtient la valeur de la variable de même indice que le premier élément non nul rencontré sur la ligne,
- les valeurs arbitraires correspondent aux zéros de la diagonale de B et en général il y a des valeurs non nulles au dessus de ces zéros, c’est pourquoi il faut initialisé la solution au vecteur des variables.

Le programme

Voici le programme de résolution d’un système linéaire :

//Veq vecteur des equations 
//v vecteur de variables
//renvoie le vecteur solution
linsolv(Veq,v):={
  local A,B,j,k,l,deq,d,res,ll,rep;
  d:=size(v);
  deq:=size(Veq);
  //A est la matrice du systeme +le 2nd membre
  A:=syst2mat(Veq,v);
  //B matrice reduite de Gauss-jordan
  B:=rref(A);
  res:=v;
  //ll ligne l de B
  ll:=makelist(0,0,d);
  for (l:=0; l<deq;l++){
    for (k:=0;k<d+1;k++){
      ll[k]:=B[l][k];
    }
  j:=l;
  while (ll[j]==0 && j<d){
    j:=j+1;
  }
  //si (j==d and ll[d]==0) 
  //ll=ligne de zeros on ne fait rien
  if (j==d and ll[d]!=0){
    // pas de sol
    return [];
  } 
  else {//la sol res[j] vaut rep/ll[j]
      if (j<d) {
      rep:=-ll[d];
      for (k:=j+1;k<d;k++) {
      rep:=rep-ll[k]*v[k];
      }
      res[j]:=rep/ll[j];
      }
    }
  }
  return res;
}

Autre algorithme

14.2.2  Résolution de MX=b donné sous forme matricielle

L’algorithme

On transforme la résolution de MX=b en MX−b=AY=0 où A est la matrice constituèe de la matrice M du système augmentèe d’une colonne formée par l’opposé du second membre b.
Y est donc un vecteur du noyau de A ayant comme dernière composante 1.
D’après l’algorithme de recherche du noyau, seul le dernier vecteur a comme dernière composante -1. Si -ker(A) renvoie n vecteurs, la solution Y est donc une combinaison arbitraire des n−1 premiers vecteurs du noyau plus le n-ième vecteur.

Le programme

Voici le programme de résolution de AX=b :

//M*res=b res renvoie le vecteur solution
//M:=[[1,1],[1,-1]]; b:=[4,2]
//M:=[[1,1,1],[1,1,1],[1,1,1]];b:=[0,0,0]
//M:=[[1,2,1],[1,2,5],[1,2,1]];b:=[0,1,0]
linsolvm(M,b):={
  local A,B,N,n,k,d,res;
  d:=ncols(M);
  //A est la matrice du systeme +le 2nd membre
  A:=border(M,-b);
  //N contient une base du noyau
  N:=-ker(A);
  n:=size(N);
  //res a d+1 composante (la derniere=1)
  res:=makelist(0,0,d);
  //C_(k) designe les constantes arbitraires 
  for (k:=0;k<n-1;k++){
    res:=res+N[k]*C_(k);
  }
  res:=res+N[n-1];
  res:=suppress(res,d);
  return res;
}

14.3  La décomposition LU d’une matrice

C’est l’interprétation matricielle de la méthode de Gauss.
Si A est une matrice carrée d’ordre n, il existe une matrice L triangulaire supérieure, une matrice L triangulaire inférieure, et une matrice P de permutation telles que :
P*A=L*U.
Supposons tout d’abord que l’on peut faire la méthode de Gauss sans échanger des lignes. Mettre des zéros sous la diagonale de la 1-ière colonne (d’ndice 0) de A, reviens à multiplier A par la matrice E₀ qui a des 1 sur sa diagonale, comme première colonne :
[1,−A[1,0]/A[0,0], ...−A[n−1,0]/A[0,0]] et des zéros ailleurs.
Puis si A1=E1*A, mettre des zéros sous la diagonale de la 2-ière colonne (d’ndice 1) de A1, reviens à multiplier A1 par la matrice E₁ qui a des 1 sur sa diagonale, comme deuxième colonne :
[0,1,−A1[2,1]/A[1,1], ...−A[n−1,1]/A[1,1]] et des zéros ailleurs.
On continue ainsi jusqu’à mettre des zéros sous la diagonale de la colonne d’indice n−2, et à la fin la matrice U=E_n−2*...*E₁*E₀*A est triangulaire supérieure et on a L=inv(E_n−2*...*E₁*E₀).
Le calcul de inv(L) est simple car on a :
- inv(E₀) des 1 sur sa diagonale, comme colonne d’indice 0 [1,+A[1,0]/A[0,0], ...+A[n−1,0]/A[0,0]] et des zéros ailleurs de même inv(E_k) est obtenue à partir de E_k "en changeant les moins en plus".
- la colonne d’indice k de inv(E₀)*inv(E₁)...inv(E_n−2) est égale à la colonne d’indice k de E_k.
Lorsqu’il y a à faire une permutation de lignes, il faut répercuter cette permutation sur L et sur U : pour faciliter la programmation on va conserver les valeurs de L et de U dans une seule matrice R que l’on séparera à la fin : U sera la partie supérieure et la diagonale de R L sera la partie inférieure de R, avec des 1 sur sa diagonale.
Voici le programme de séparation :

splitmat(R):={
  local L,U,n,k,j;
  n:=size(R);
  L:=idn(n);
  U:=makemat(0,n,n);
  for (k:=0;k<n;k++){
    for (j:=k;j<n;j++){
      U[k,j]:= R[k,j];
    }
  } 
  for (k:=1;k<n;k++){
    for (j:=0;j<k;j++){
      L[k,j]:= R[k,j];
    }
  }
return (L,U);
};

Le programme ci-dessous, decomplu(A), renvoie la permutation p que l’on a fait sur les lignes, L et U et on a: P*A=L*U.
Voici le programme de décomposition LU qui utilise splitmat ci-dessus :

//A:=[[5,2,1],[5,2,2],[-4,2,1]]
//A:=[[5,2,1],[5,-6,2],[-4,2,1]]
// utilise splitmat
decomplu(A):={
  local B,R,L,U,n,j,k,l,temp,p;
  n:=size(A);
  p:=seq(k,k,0,n-1);
  R:=makemat(0,n,n);
  B:=A;
  l:=0;
//on traite toutes les colonnes
  while (l<n-1) {
    if (A[l,l]!=0){
    //pas de permutations 
    //on recopie dans R la ligne l de A 
    //a partir de la diagonale
      for (j:=l;j<n;j++){R[l,j]:=A[l,j];}
      //on met des zeros sous la diagonale 
      //dans la colonne l
      for (k:=l+1;k<n;k++){
 for (j:=l+1;j<n;j++){
   A[k,j]:=A[k,j]-A[l,j]*A[k,l]/A[l,l];
   R[k,j]:=A[k,j];
 } 
 R[k,l]:=A[k,l]/A[l,l];
 A[k,l]:=0;
      }
    l:=l+1;
    }
    else {
      k:=l;
      while ((k<n-1) and (A[k,l]==0)) {
        k:=k+1;
      }
      //si (A[k,l]==0) A est non inversible, 
      //on passe a la colonne suivante
      if (A[k,l]==0) {
 l:=l+1;
      }
      else {
      //A[k,l]!=0) on echange la ligne l et k ds A et R
       for (j:=l;j<n;j++){
        temp:=A[k,j];
        A[k,j]:=A[l,j];
        A[l,j]:=temp;
       };
       for (j:=0;j<n;j++){
        temp:=R[k,j];
        R[k,j]:=R[l,j];
        R[l,j]:=temp;
       }
       //on note cet echange dans  p
       temp:=p[k];
       p[k]:=p[l];
       p[l]:=temp;
     }//fin du if (A[k,l]==0)
    }//fin du if (A[l,l]!=0)
  }//fin du while
  L,U:=splitmat(R);
  return(p,L,U);
};

On tape :
A:=[[5,2,1],[5,2,2],[-4,2,1]]
decomplu(A)
On obtient :
[0,2,1],[[1,0,0],[-4/5,1,0],[1,0,1]],
[[5,2,1],[0,18/5,9/5],[0,0,1]]
On verifie, on tape (car permu2mat(p)=P=inv(P)) :
[[1,0,0],[0,0,1],[0,1,0]]*[[1,0,0],[-4/5,1,0],[1,0,1]]*
[[5,2,1],[0,18/5,9/5],[0,0,1]]
On obtient bien [[5,2,1],[5,2,2],[-4,2,1]]
On tape :
B:=[[5,2,1],[5,-6,2],[-4,2,1]]
decomplu(B)
On obtient :
[0,1,2],,[[1,0,0],[1,1,0],[-4/5,-9/20,1]],
[[5,2,1],[0,-8,1],[0,0,9/4]]
On verifie, on tape :
[[1,0,0],[1,1,0],[-4/5,-9/20,1]]*[[5,2,1],[0,-8,1],[0,0,9/4]]
On obtient bien [[5,2,1],[5,-6,2],[-4,2,1]]

14.4  Décomposition de Cholesky d’une matrice symétrique définie positive

14.4.1  Les méthodes

Lorsque la matrice A est la matrice associée à une forme bilinéaire définie positive, on lui associe la matrice symétrique B=1/2*(A+tran(A)) qui est la matrice de la forme quatratique associée.
Avec ses hypothèses, il existe une matrice triangulaire inférieure unique C telle que B=C*tran(C).
Pour déterminer C on présente ici deux méthodes :

La méthode utilisant la décomposition LU

Si on décompose B selon la méthode LU, on n’est pas obligé de faire des échanges de lignes car les sous-matrices principales B_k d’ordre k (obtenues en prenant les k premières lignes et colonnes de B) sont des matrices inversibles car ce sont des matrices de formes définies positives.
On a donc :
p,L,U:=decomplu(B) avec p=[0,1..n-1] si A est d’ordre n. Posons D la matrice diagonale ayant comme diagonale la racine carrée de la diagonale de U.
On a alors C=L*D et tran(C)=inv(D)*U.
Pour le montrer on utilise le théorème :
Si B=C*F avec B symétrique, C triangulaire inférieure et F triangulaire supérieure de même diagonale que C alors F=tran(C).
En effet on a :
B=tran(B) donc C*F=tran(C*F)=tran(F)*tran(C).
On en déduit l’égalité des 2 matrices :
inv(tran(F))*C=tran(C)*inv(F)
la première est triangulaire inférieure, et la deuxième est triangulaire supérieure donc ces matrices sont diagonales et leur diagonale n’a que des 1 ces 2 matrices sont donc égales à la matrice unité.

La méthode par identification

On peut aussi déterminer C par identification en utilisant l’égalité B=C*tran(C) et C triangulaire inférieure.
On trouve pour tout entier j compris entre 0 et n−1 :
- pour la diagonale :
(C[j,j])²=B[j,j]-∑_k=0^j-1(C[j,k])²
- pour les termes subdiagonaux c’est à dire pour j+1 ≤ l ≤ n−1 :
C[l,j]=1/C[j,j]*(B[l,j]-∑_k=0^j-1C[l,k]*C[j,k])
On peut donc calculer les coefficients de C par colonnes, en effet, pour avoir la colonne j :
- on calcule le terme diagonal C[j,j] qui fait intervenir les termes de la ligne j des colonnes précédentes (avec un test vérifiant que le terme B[j,j]-∑_k=0^j-1(C[j,k])² est positif), puis,
- on calcule les termes subdiagonaux C[l,j] pour j+1 ≤ l ≤ n−1 qui font intervenir les termes des colonnes précédentes des lignes précédentes. Mais cet algorithme n’est pas très bon car les termes de la matrice C sont obtenus à chaque étape avec une multiplication ou une division de racine carrée : il serait préférable d’introduire les racines carrées qu’à la fin du calcul comme dans la méthode LU. On va donc calculer une matrice CC (sans utiliser de racines carrées) et une matrice diagonale D (qui aura des racines carrées sur sa diagonale) de tel sorte que C=CC*D (la colonne k de CC*D est la colonne k de CC multiplié par D[k,k]).
Par exemple : C[0,0]* C[0,0]= B[0,0] on pose :
CC[0,0]= B[0,0]=a et D[0,0]=1/sqrt(a) ainsi
C[0,0]= CC[0,0]/sqrt(a)=sqrt(a) et on a bien :
C[0,0]* C[0,0]=a=B[0,0]
On aura donc C[l,0]=CC[l,0]/sqrt(a) c’est à dire :
CC[l,0]=B[l,0] pour 0≤ l<n
Puis :
C[1,1]*C[1,1]=B[1,1]-C[1,0]*C[1,0]=
B[1,1]-CC[1,0]*CC[1,0]/a=b
on pose :
CC[1,1]= B[1,1]-CC[1,0]*CC[1,0]/a=b et
D[1,1]=1/sqrt(b)
On aura donc pour 2≤ l<n ;
C[l,1]=CC[l,1]/sqrt(b)=1/sqrt(b)*(B[l,1]-C[l,0]*C[1,0])
c’est à dire pour 2≤ l<n :
CC[l,1]=B[l,1]-C[l,0]*C[1,0]=B[l,1]-CC[l,0]*CC[1,0]/a
etc.....
Les formules de récurences pour le calcul de CC sont :
- pour la diagonale :
CC[j,j]=B[j,j]-∑_k=0^j-1(CC[j,k])²/CC[k,k]
- pour les termes subdiagonaux c’est à dire pour j+1 ≤ l ≤ n−1 :
CC[l,j]=B[l,j]-∑_k=0^j-1CC[l,k]*CC[j,k]/CC[k,k]
avec D[j,j]=1/sqrt(CC[j,j]).

14.4.2  Le programme de factorisation de Cholesky avec LU

 
//utilise decomplu et splitmat ci-dessus
//A:=[[1,0,-1],[0,2,4],[-1,4,11]]
//A:=[[1,1,1],[1,2,4],[1,4,11]]
//A:=[[1,0,-2],[0,2,6],[0,2,11]]
//A:=[[1,-2,4],[-2,13,-11],[4,-11,21]]
//A:=[[-1,-2,4],[-2,13,-11],[4,-11,21]] (pas def pos)
//A:=[[24,66,13],[66,230,-11],[13,-11,210]]
choles(A):={
  local j,n,p,L,U,D,p0;
  n:=size(A);
  A:=1/2*(A+tran(A));
 (p,L,U):=decomplu(A);
  p0:=makelist(x->x,0,n-1);
  if (p!=p0) {return "pas definie positive ";}
  D:=makemat(0,n,n);
  for (j:=0;j<n;j++) {
    //if (U[j,j]<0) {return "pas def positive";}
    D[j,j]:=sqrt(U[j,j]);
  }
  return normal(L*D);
}

14.4.3  Le programme de factorisation de Cholesky par identification

 
//A:=[[1,0,-1],[0,2,4],[-1,4,11]]
//A:=[[1,1,1],[1,2,4],[1,4,11]]
//A:=[[1,0,-2],[0,2,6],[0,2,11]]
//A:=[[1,-2,4],[-2,13,-11],[4,-11,21]]
//A:=[[-1,-2,4],[-2,13,-11],[4,-11,21]] (pas def pos)
//A:=[[24,66,13],[66,230,-11],[13,-11,210]]
cholesi(A):={
  local j,n,l,k,C,c2,s;
  n:=size(A);
  A:=1/2*(A+tran(A));
  C:=makemat(0,n,n);;
  for (j:=0;j<n;j++) {
    s:=0;
    for (k:=0;k<j;k++) {
      s:=s+(C[j,k])^2;
    }
    c2:=A[j,j]-s; 
    if (c2<=0) {return "pas definie positive ";}
    C[j,j]:=normal(sqrt(c2));
    for (l:=j+1;l<n;l++) {
      s:=0;
      for (k:=0;k<j;k++) {
 s:=s+C[l,k]*C[j,k];
      }
      C[l,j]:=normal(1/sqrt(c2)*(A[l,j]-s));
    }
  }
return C;
}

Ce programme n’est pas très bon car les termes de la matrice C sont obtenus avec une multiplication ou une division de racine carrée...
On se pourra se reporter ci-dessous pour avoir le programme choleski optimisé.

14.4.4  Le programme optimisé de factorisation de Cholesky par identification

 
//A:=[[1,0,-1],[0,2,4],[-1,4,11]]
//A:=[[1,1,1],[1,2,4],[1,4,11]]
//A:=[[1,0,-2],[0,2,6],[0,2,11]]
//A:=[[1,-2,4],[-2,13,-11],[4,-11,21]]
//A:=[[-1,-2,4],[-2,13,-11],[4,-11,21]] (pas def pos)
//A:=[[24,66,13],[66,230,-11],[13,-11,210]]
choleski(A):={
  local j,n,l,k,CC,D,c2,s;
  n:=size(A);
  A:=1/2*(A+tran(A));
  CC:=makemat(0,n,n);
  D:=makemat(0,n,n);
  for (j:=0;j<n;j++) {
    s:=0;
    for (k:=0;k<j;k++) {
      s:=s+(CC[j,k])^2/CC[k,k];
    }
    c2:=normal(A[j,j]-s); 
    //if (c2<=0) {return "pas definie positive ";}
    CC[j,j]:=c2;
    D[j,j]:=normal(1/sqrt(c2));
    for (l:=j+1;l<n;l++) {
      s:=0;
      for (k:=0;k<j;k++) {
 s:=s+CC[l,k]*CC[j,k]/CC[k,k];
      }
      CC[l,j]:=normal(A[l,j]-s);
    }
  }
return normal(CC*D);
}

Avec cette méthode, pour obtenir les coefficients diagonaux ont utilise la même relation de récurrence que pour les autres coefficients. De plus on peut effectuer directement et facilement la multiplication par D sans avoir à définir D en multipliant les colonnes de CC par la même valeur 1/sqrt(c) avec c=CC[k,k]...
donc on écrit :

 
choleskii(A):={
  local j,n,l,k,CC,c,s;
  n:=size(A);
  A:=1/2*(A+tran(A));
  CC:=makemat(0,n,n);
  for (j:=0;j<n;j:=j+1) { 
     for (l:=j;l<n;l++) {
      s:=0;
      for (k:=0;k<j;k++) {
 //if (CC[k,k]<=0) {return "pas definie positive ";}
 if (CC[k,k]==0) {return "pas definie";}
 s:=s+CC[l,k]*CC[j,k]/CC[k,k];
      }
      CC[l,j]:=A[l,j]-s;
    }
  } 
  for (k:=0;k<n;k++) {
    c:=CC[k,k];
    for (j:=k;j<n;j++) {
      CC[j,k]:=normal(CC[j,k]/sqrt(c));
    }
  }
  return CC;
}

Avec la traduction pour avoir une fonction interne à Xcas :

 
cholesky(_args)={
gen args=(_args+mtran(_args))/2;
matrice &A=*args._VECTptr;
int n=A.size(),j,k,l;
vector<vecteur> C(n,vecteur(n)), D(n,vecteur(n));
for (j=0;j<n;j++) {
gen s;
for (k=0;k<j;k++)
s=s+pow(C[j][k],2)/C[k][k];
gen c2=A[j][j]-s;
if (is_strictly_positive(-c2)) 
setsizeerr("Not a positive define matrice");
C[j][j]=c2;
D[j][j]=normal(1/sqrt(c2));
for (l=j+1;l<n;l++) {
s=0;
for (k=0;k<j;k++) 
s=s+C[l][k]*C[j][k]/C[k][k];
C[l][j]=A[l][j]-s;
}
}
matrice Cmat,Dmat;
vector_of_vecteur2matrice(C,Cmat);
vector_of_vecteur2matrice(D,Dmat);
return Cmat*Dmat;
}

14.5  Réduction de Hessenberg

14.5.1  La méthode

Une matrice de Hessenberg est une matrice qui a des zéros sous la "deuxième diagonale inférieure".
Soit A une matrice. On va chercher B une matrice de Hessenberg semblable à A. Pour cela on va mettre des zéros sous cette diagonale en utilisant la méthode de Gauss mais en prenant les pivots sur la "deuxième diagonale inférieure" encore appelée "sous-diagonale": cela revient à multiplier A par Q=R⁻¹ et cela permet de conserver les zéros lorsque l’on multiplie à chaque étape le résultat par R pour obtenir une matrice semblable à A. Si on est obligé de faire un échange de lignes (correspondant à la multiplication à droite par E) il faudra faire aussi un échange de colonnes (correspondant à la multiplication à gauche par E.
Par exemple si on a :
A:=[

]
Pour mettre des zéros dans la première colonne en dessous de a₁₀, on va multiplier A à gauche par Q et à droite par R=Q⁻¹ avec si on suppose a₁₀!=0 c’est à dire si on peut choisir comme pivot a₁₀ : Q:=[

]
R:=[

]
On tape alors :
A:=[[a₀₀,a₀₁,a₀₂,a₀₃,a₀₄],[a₁₀,a₁₁,a₁₂,a₁₃,a₁₄],
[a₂₀,a₂₁,a₂₂,a₂₃,a₂₄],[a₃₀,a₃₁,a₃₂,a₃₃,a₃₄],[a₄₀,a₄₁,a₄₂,a₄₃,a₄₄]]
Q:=[[1,0,0,0,0],[0,1,0,0,0],[0,(-a₂₀)/a₁₀,1,0,0],
[0,(-a₃₀)/a₁₀,0,1,0],[0,(-a₄₀)/a₁₀,0,0,1]]
R:=[[1,0,0,0,0],[0,1,0,0,0],[0,(a₂₀)/a₁₀,1,0,0],
[0,(a₃₀)/a₁₀,0,1,0],[0,(a₄₀)/a₁₀,0,0,1]]
On obtient la matrice B1:
B1=Q*A*R=R⁻¹*A*R=[

] où les ... sont des expressions des coefficients de A.
On va faire cette transformation successivement sur la deuxième colonne de la matrice B1 pour obtenir B2 etc...
On appellera B toutes les matrices obtenues successivement.
Si on doit échanger les deux lignes k et j pour avoir un pivot, cela revient à multiplier à gauche la matrice B par une matrice E égale à la matrice identité ayant subie l’échange des deux lignes k et j. Il faudra alors, aussi multiplier à droite la matrice B par E c’est à dire échanger les deux colonnes k et j de B

14.5.2  Le programme de réduction de Hessenberg

 
//A est mise sous forme de hessenberg B avec 
//P matrice de passage
//p indique si il y a eu une permutation de lignes
//a chaque etape la matrice inv(R) 
// met des zeros sous la "sous diagonale"
//B=P^-1AP ex A:=[[5,2,1],[5,2,2],[-4,2,1]]
//A:=[[5,2,1,1],[5,2,2,1],[5,2,2,1],[-4,2,1,1]]
hessenbg(A):={
  local B,R,P,n,j,k,l,temp,p;
  n:=size(A);
  P:=idn(n);
  p:=seq(k,k,0,n-1);
  B:=A;
  l:=1;
  while (l<n) {
     R:=idn(n);
     if (B[l,l-1]!=0){
     for (k:=l+1;k<n;k++){
 R[k,l]:=B[k,l-1]/B[l,l-1];
 //on multiplie B a droite par inv(R)
        for (j:=l;j<n;j++){
   B[k,j]:=B[k,j]-B[l,j]*B[k,l-1]/B[l,l-1];
 } 
 B[k,l-1]:=0;
      }
     //on multiplie B et P a gauche par R
     B:=B*R;
     P:=P*R;
    l:=l+1;
    }
  else {
    k:=l;
    while ((k<n-1) and (B[k,l-1]==0)) {
    k:=k+1;
    }
    
    if (B[k,l-1]==0) {l:=l+1;}
    else{
    //B[k,l]!=0) on echange ligne l et k ds B
    for (j:=l-1;j<n;j++){
    temp:=B[k,j];
    B[k,j]:=B[l,j];
    B[l,j]:=temp;
    };
    //A[k,l]!=0) on echange colonne l et k ds B et P
    for (j:=0;j<n;j++){
    temp:=B[j,k];
    B[j,k]:=B[j,l];
    B[j,l]:=temp;
    };
    for (j:=0;j<n;j++){
    temp:=P[j,k];
    P[j,k]:=P[j,l];
    P[j,l]:=temp;
    };
    temp:=p[k];
    p[k]:=p[l];
    p[l]:=temp;
  }
  }
  }
  return(p,P,B);
};

p nous dit les échanges effectués et on a B=P⁻¹*A*P.
Attention !! si A est symétrique, cette transformation détruit la symétrie et on n’a donc pas une tridiagonalistion d’une matrice symtrisue par cette méthode.

14.6  Tridiagonalisation des matrices symétriques avec des rotations

On a le théorème : Pour toute matrice symétrique A d’ordre n, il existe une matrice P, produit de n−2 matrices de rotations, telle que B=^tP*A*P soit tridiagonale : c’est la réduction de Givens.

14.6.1  Matrice de rotation associée à e_p, e_q

Dans ℝⁿ, on appelle rotation associée à e_p, e_q, une rotation d’angle t du plan dirigé par e_p, e_q où e_k désigne le k+1-ième vecteur de la base canonique de ℝⁿ (la base canonique est e₀=[1,0..,0], e₁=[0,1,0..,0] etc...).
Si ℝⁿ est rapporté à la base canonique, à cette rotation est associée une matrice G(n,p,q,t) dont le terme général est :
si k∉{p,q}, g_k,k=1
g_p,p=g_q,q=cos(t)
g_p,q=−g_q,p=−sin(t)
Voici un matrice de rotation associée à e₁, e₃ (deuxième et quatrième vecteur de base) de ℝ⁵ :

et le programme qui construit une telle matrice :

rota(n,p,q,t):={
local G;
G:=idn(n);
G[p,p]:=cos(t);
G[q,q]:=cos(t);
G[p,q]:=-sin(t);
G[q,p]:=sin(t);
return G;
}

14.6.2  Réduction de Givens

Soit A une matrice symétrique. On va chercher G1 une matrice de rotation associée à e₁,e₂ pour annuler les coefficients 2,0 et 0,2 de ^tG1*A*G1.
Regardons un exemple :
G1:=[[1,0,0,0,0],[0,cos(t),-sin(t),0,0],[0,sin(t),cos(t),0,0],
[0,0,0,1,0],[0,0,0,0,1]]
A:=[[a,b,c,d,e],[b,f,g,h,j],[c,g,k,l,m],[d,h,l,n,o],[e,j,m,o,r]]
On obtient :
^tG1=[

], A= [

]
On a :
^tG1*A= [

]
On choisit t pour que −bsin(t)+ccos(t)=0 par exemple :
cos(t)=b/√(b²+c²) et sin(t)=c/√(b²+c²).
On a :
G1:=[

]
et donc puisqu’on a choisit −bsin(t)+ccos(t)=0, on a bien annuler les coefficients 2,0 et 0,2 de :
^tG1*A*G1= [

]
Puis on annule les coefficients 0,3 et 3,0 de l la matrice A1 obtenue en formant ^tG2*A1*G2 avec :
G2:=[

] en choissant correctement t etc...

14.6.3  Le programme de tridiagonalisation par la méthode de Givens

//A:=[[1,-1,2,1],[-1,1,-1,1],[2,-1,1,-1],[1,1,-1,-1]]
//tran(R)*A*R=B si tridiagivens(A)=R,B
//pour annuler le terme 2,0 on multiplie A par TG
//TG[1,1]=cos(t)=TG[2,2],TG[1,2]=sin(t)=-TG[2,1], 
//avec -sin(t)A[1,0]+cos(t)A[2,0]=0
//TG*A multiplie par tran(TG) annule les termes 2,0 et 0,2
//pour annuler le terme q,p (q>p+1) on multiplie A par TG
//TG[p+1,p+1]=cos(t)=TG[q,q],TG[p+1,q]=sin(t)=-TG[q,p+1], 
//avec sin(t)A[p+1,p]=cos(t)A[q,p] 
//donc sin(t)=A[q,p]/r et cos(t)=A[p+1,p]/r
//avec r:=sqrt((A[p+1,p])^2+(A[q,p])^2)
tridiagivens(A):={
  local n,p,q,r,TG,R,c,s;
  n:=size(A);
  R:=idn(n);
  for (p:=0;p<n-2;p++) {
    for (q:=p+2;q<n;q++) {
    r:=sqrt((A[p+1,p])^2+(A[q,p])^2);
    if (r!=0) {
    c:=normal(A[p+1,p]/r);
    s:=normal(A[q,p]/r);
    TG:=idn(n);
    TG[p+1,p+1]:=c;
    TG[q,q]:=c;
    TG[p+1,q]:=s;    
    TG[q,p+1]:=-s;
    TG:=normal(TG);
    A:=normal(TG*A);
    A:=normal(A*tran(TG));
    A:=normal(A);
    R:=normal(R*tran(TG));
    }
    }
  }
  return (R,A);
}   

On tape :
A:=[[1,-1,2,1],[-1,1,-1,1],[2,-1,1,-1],[1,1,-1,-1]]
tridiagivens(A)
On obtient :

[[1,0,0,0],[0,(-(sqrt(6)))/6,
(-(sqrt(4838400000)))/201600,
(-(sqrt(2962400000)))/64400],
[0,sqrt(6)/3,sqrt(4838400000)/252000,
(-(sqrt(26661600000)))/322000],
[0,sqrt(6)/6,(-(26*sqrt(210)))/420,12*sqrt(35)/420]],
[[1,sqrt(6),0,0],
[sqrt(6),1/3,sqrt(275990400000)/266400,0],
[0,sqrt(209026944000)/231840,73/105,8*sqrt(6)/35],
[0,0,8*sqrt(6)/35,1/-35]]

On tape :
A:=[[1,-1,2,0],[-1,1,-1,1],[2,-1,1,-1],[1,1,-1,-1]]
tridiagivens(A)
On obtient :

[[1,0,0,0],[0,(-(sqrt(6)))/6,
(-(sqrt(4838400000)))/201600,(-(sqrt(2962400000)))/64400],
[0,sqrt(6)/3,sqrt(4838400000)/252000,
(-(sqrt(26661600000)))/322000],
[0,sqrt(6)/6,(-(26*sqrt(210)))/420,12*sqrt(35)/420]],
[[1,5*sqrt(6)/6,13*sqrt(210)/210,(-(sqrt(35)))/35],
[sqrt(6),1/3,sqrt(275990400000)/266400,0],
[0,sqrt(209026944000)/231840,73/105,8*sqrt(6)/35],
[0,0,8*sqrt(6)/35,1/-35]]

14.7  Tridiagonalisation des matrices symétriques avec Householder

On a le théorème : Pour toute matrice symétrique A d’ordre n, il existe une matrice P, produit de n−2 matrice de Householder (donc P⁻¹=^tP), telle que B=^tP*A*P soit tridiagonale.
Si A n’est pas symétrique, il existe une matrice P, produit de n−2 matrice de Householder, telle que B= ^tP*A*P soit un matrice de Hessenberg c’est à dire une matrice dont les coefficient sous-tridiagonaux sonts nuls : c’est la réduction de Householder.
On va dans cette section écrire un programme qui renverra P et B.

14.7.1  Matrice de Householder associée à v

Soit un hyperplan défini par son vecteur normal v.
On appelle matrice de Householder, la matrice d’une symétrie orthogonale hv par rapport à cet hyperplan. On a :
hv(u):=u-2(dot(v,u))/norm(v)^2*v
et la matrice de Householder associée est :
Hv=idn(n)-2*tran(v)*[v]/(v*v).
On a :
Hv=tran(Hv)=inv(Hv) Par convention la matrice unité est une matrice de Householder.
On écrit le programme houdeholder qui à un vecteur v renvoie la matrice de Householder associé à v.

householder(v):={
  local n,w,nv;
  n:=size(v);
  nv:=v*v;
  w:=tran(v);
  return normal(idn(n)-2*w*[v]/nv);
};

14.7.2  Matrice de Householder annulant les dernières composantes de a

Soit Hv la matrice de Householder associée à v.
Etant donné un vecteur a non nul, il existe deux vecteurs v tels que le vecteur b=Hv*tran(a) a ses n−1 dernières composantes nulles (Hv étant la matrice de Householder associé à v).
Plus précisement, si e0=[1,0,..0],e1=[0,1,0,..0] etc..., les vecteurs v sont :
v1=a+norm(a)*e0 et v2=a-norm(a)*e0.
Plus généralement, il existe deux vecteurs v tels que le vecteur b=Hv*tran(a) a ses composantes à partir de la (l+1)-ième nulles.
Plus précisement, si el=[0,..0,1,0,..0] (ou toutes les composantes sont nulles sauf el[l]=1), et si al=[0,...,0,a[l],..,a[n-1]], les vecteurs v sont :
v1=al+norm(al)*el et v2=al-norm(al)*el.
On 'ecrit le programme qui étant donné a et l renvoie si a est nul renvoie la matrice identité et si a est non nul renvoie la matrice de Householder associée à v=al+eps*norm(al)*el en choisisant eps égal au signe de al si al est réelle et eps=1 sinon.

//renvoie une matrice hh= H tel que b:=H*tran(a) 
//a ses n-l-1 dernieres comp=0 
//ou encore b[l+1]=...b[n-1]=0
make_house(a,l):={
  local n,na,el,v;
  n:=size(a);
  for (k:=0;k<l;k++) a[k]:=0;
  na:=normal(norm(a)^2);
  na:=sqrt(na);
  if (na==0) return idn(n);
  el:=makelist(0,1,n);
  el[l]:=1;
  if (im(a[l])==0 and a[l]<0) v:=normal(a-na*el); 
    else v:=normal(a+na*el);
  return householder(v);
};

14.7.3  Réduction de Householder

Le programme tridiaghouse(A) renvoie H,B tel que normal(H*A*tran(H))=B (ou encore normal(tran(H)*B*H)=A) avec H est le produit de matrice de householder (tran(H)=inv(H)) et B est une matrice de hessenberg (ou bien B est une matrice tridiagonale si tran(A)=A).

tridiaghouse(A):={
  local n, a, H,Ha;
  B:=A;
  n:=size(A);
  H:=idn(n);
  for (l:=0;l<n-2;l++) {
  a:=col(B,l);
  Ha:=make_house(a,l+1);
  B:=normal(Ha*B*Ha);
  H:=normal(Ha*H);
  }
  //normal(H*A*tran(H))=B et B de hessenberg 
  //(ou tridiagonale si tran(A)=A)
  return(H,B);
};

On tape :
A:=[[3,2,2,2,2],[2,1,2,-1,-1],[2,2,1,-1,1],
[2,-1,-1,3,1],[2,-1,1,1,2]]
H,B:tridiaghouse(A)
On obtient :
[[1,0,0,0,0],[0,1/-2,1/-2,1/-2,1/-2],[0,5*sqrt(11)/22,
(-(3*sqrt(11)))/22,sqrt(11)/22,(-(3*sqrt(11)))/22],
[0,0,22*sqrt(2)/44,0,(-(22*sqrt(2)))/44],[0,22*sqrt(22)/242,
22*sqrt(22)/484,(-(22*sqrt(22)))/121,22*sqrt(22)/484]],
[[3,-4,0,0,0],[-4,9/4,sqrt(11)/4,0,0],[0,sqrt(11)/4,3/4,
44*sqrt(22)/121,0],[0,0,44*sqrt(22)/121,1/2,286*sqrt(11)/484],
[0,0,0,286*sqrt(11)/484,7/2]]
On vérifie et on tape :
normal(H*A*tran(H))
On tape :
A:=[[1,2,3],[2,3,4],[3,4,5]]
H,B:=tridiaghouse(A)
On obtient :
[[1,0,0],[0,(-(2*sqrt(13)))/13,(-(3*sqrt(13)))/13],
[0,(-(3*sqrt(13)))/13,(2*sqrt(13))/13]],
[[1,-(sqrt(13)),0],[-(sqrt(13)),105/13,8/13],[0,8/13,(-1)/13]]
On vérifie et on tape :
normal(H*A*tran(H))