Inversion de matrice par la méthode de Pan et Reiff.


Un peu de théorie

Matrice

On appelle matrice un tableau rectangulaire de nombres réels ou complexes. La méthode d'inversion décrite ci-dessous le sera dans le cadre des nombres réels mais il en existe une variante pour les complexes.

Exemple:
    -1,2  2,5     0
    12,1  0,1  -2,3
La matrice ci-dessus possède 6 éléments rangés dans 2 lignes de 3 colonnes. Le symbole utilisé pour la représenter est A ou plus précisément A2x3 où 2x3 indique qu'elle possède 2 lignes et 3 colonnes.
Les éléments d'une matrice sont représentés par une minuscule (a pour la matrice A) accompagnée d'indices pour indiquer la position de l'élément: ai,j où i représente le numéro de la ligne et j le numéro de la colonne.
Dans la matrice de l'exemple, a2,3 vaut -2,3 et a1,2 vaut 2,5.
Seules les matrices carrées, même nombre de lignes et de colonnes, peuvent être éventuellement inversées.

Somme de matrices

Seules sont permises les sommes/différences de matrices de même dimension.

On aura ANxN + BNxN = CNxN avec
    ci,j = ai,j + bi,j pour i et j=1, ..., N

Produit de matrices

Seuls sont permis les produits AxB de matrices compatibles par leur dimension.

En effet, l'élément ci,j du produit étant la somme des éléments de la ième ligne de A et de la jème colonne de B, il faut que les lignes de A et les colonnes de B possèdent le même nombre d'éléments.

On aura: ALxM + BMxN = CLxN avec
    ci,j = Somme (ai,k*bk,j) pour k=1, .., M

Matrice inverse

On appelle inverse de la matrice A la matrice notée A-1 telle que
    A-1 . A = A . A-1 = I
où la matrice I est la matrice identité de même dimension que A et composée de 0 sauf les éléments diagonaux qui valent 1.

La méthode utilisée pour inverser une matrice est souvent:

* transposer (voir plus bas) la matrice adjointe (matrice des cofacteurs transposée)
* diviser chaque élément par le déterminant

ce qui revient à:
1] calculer la comatrice (la matrice des cofacteurs)
2] diviser chaque élément par le déterminant.

Ceci montre que la condition nécessaire et suffisante pour qu'une matrice carrée soit inversible est que son déterminant soit non nul.

Le réel problème est donc dans le calcul du déterminant.
* La méthode des mineurs/cofacteurs déjà peu pratique pour plus de 3 lignes (répétition de mêmes calculs) quand il s'agit de l'utiliser "à la main" n'est pas plus facile à programmer.
* La méthode d'élimination de Gauss, plus simple à programmer, nécessite malheureusement de nombreux calculs (de l'ordre de n3) et est donc peu utilisée pour n grand.

Matrice transposée

On appelle matrice transposée d'une matrice A et on la note At la matrice obtenue en prenant l'élément de A en ligne i, colonne j et en le plaçant en ligne j et colonne i.

On a donc:
    ati,j = aj,i pour i,j=1, .., N

Norme cartésienne

La norme a pour but d'essayer de quantifier par un seul nombre la taille des éléments de la matrice. Il en existe plusieurs dont la norme cartésienne qui est la racine carrée de la somme des carrés de tous les éléments de la matrice. Cette norme a le gros avantage de ne pas faire intervenir de tests ("gros" consommateurs de ressources). Cette norme peut être vue une extension du théorème de Pythagore.
    Norme = Somme (ai,j*ai,j) pour i,j=1, .., N

Norme ligne

Pour déterminer cette norme, il faut, pour chaque ligne, calculer la somme des valeurs absolues des éléments et prendre la plus grande valeur.


Norme colonne

Pour déterminer cette norme, il faut, pour chaque colonne, calculer la somme des valeurs absolues des éléments et prendre la plus grande valeur.

L'algorithme

La méthode de Pan et Reiff est une méthode itérative qui nécessite une première approximation "valable" pour que la convergence du processus soit assurée. Tout le mérite de la découverte de cette première "valeur" revient à Pan et Reiff qui ont ainsi réussi à rendre opérationnelle une méthode déjà bien connue.

Première approximation

Pour que le processus itératif soit convergent, il suffit de prendre comme première approximation la matrice transposée et de diviser chaque élément par le produit des normes ligne et colonne.

Les itérations

* Déterminer la première approximation B comme décrit ci-dessus;
* Calculer la matrice de l'erreur E, à savoir la différence entre ce que l'on devrait avoir (I=A.A-1) et ce que l'on a (A.B).
   E = I - B.A
* Calculer la nouvelle valeur de B
   B=(I+E)*B (et non B=B.(I+E) comme indiqué auparavant. Désolé!)
* Recommencer jusqu'au moment où la norme de E est "suffisamment" petite ou quand la norme augmente. Comme au début du processus, la première approximation n'est pas toujours suffisamment précise, ce critère ne doit pas être pris en considération avant la 5ème itération.

Le programme en Turbo Pascal
PROGRAM Pan_Reiff;

USES Crt;

TYPE Matrice=ARRAY[1..10,1..10] OF REAL;

VAR MatA,MatB,MatC:Matrice;
    N:BYTE;
    OK:BOOLEAN;

PROCEDURE Intro(VAR N:BYTE;VAR Mat:Matrice);
VAR i,j:BYTE;
BEGIN
CLRSCR;
TEXTCOLOR(YELLOW);
WRITE('Nombre de lignes:');READLN(N);
FOR i:=1 TO N DO
  FOR j:=1 TO N DO
    BEGIN
    {GOTOXY(3,4);CLREOL;
    WRITE('Introduisez A(',i,',',j,'): ');
    READLN(Mat[i,j]);}
    Mat[i,j]:=RANDOM(100);{Pour tester}
    END;
END;

PROCEDURE Attente;
VAR St:STRING;
BEGIN
TEXTCOLOR(LIGHTRED);
St:='Une touche pour continuer.';
GOTOXY(40-LENGTH(St) DIV 2,25);WRITE(St);READKEY;
END;

PROCEDURE Affiche(S:STRING;N:BYTE;Mat:Matrice);
VAR i,j:BYTE;
BEGIN
CLRSCR;
TEXTCOLOR(LIGHTGREEN);
WRITELN(S);
FOR i:=1 TO N DO
  BEGIN
  FOR j:=1 TO N DO
      WRITE(Mat[i,j]:9:4);
  WRITELN;
  END;
END;

PROCEDURE Somme(N:BYTE;MatA,MatB:Matrice;VAR MatC:Matrice);
VAR i,j:BYTE;
BEGIN
FOR i:=1 TO N DO
  FOR j:=1 TO N DO
     MatC[i,j]:=MatA[i,j]+MatB[i,j];
END;

PROCEDURE Difference(N:BYTE;MatA,MatB:Matrice;VAR MatC:Matrice);
VAR i,j:BYTE;
BEGIN
FOR i:=1 TO N DO
  FOR j:=1 TO N DO
     MatC[i,j]:=MatA[i,j]-MatB[i,j];
END;

PROCEDURE Produit(N:BYTE;MatA,MatB:Matrice;VAR MatC:Matrice);
VAR i,j,k:BYTE;
BEGIN
FOR i:=1 TO N DO
  FOR j:=1 TO N DO
    BEGIN
    MatC[i,j]:=0;
    FOR k:=1 TO N DO
      MatC[i,j]:=MatC[i,j]+MatA[i,k]*MatB[k,j];
    END;
END;

PROCEDURE Inverse(N:BYTE;A:Matrice;VAR X:Matrice;VAR OK:BOOLEAN);
VAR Id,E,B:Matrice;
    Old,New:REAL;
    Cpt:BYTE;

  FUNCTION Norme_Ligne(N:BYTE;B:Matrice):REAL;
  VAR Tmp,Tmp2:REAL;
      i,j:BYTE;
  BEGIN
  Tmp:=0;
  FOR i:=1 TO N DO
    BEGIN
    Tmp2:=0;
    FOR j:=1 TO N DO
      Tmp2:=Tmp2+ABS(B[i,j]);
    IF Tmp2>Tmp THEN Tmp:=Tmp2;
    END;
  Norme_Ligne:=Tmp
  END;

  FUNCTION Norme_Colonne(N:BYTE;B:Matrice):REAL;
  VAR Tmp,Tmp2:REAL;
      i,j:BYTE;
  BEGIN
  Tmp:=0;
  FOR i:=1 TO N DO
    BEGIN
    Tmp2:=0;
    FOR j:=1 TO N DO
      Tmp2:=Tmp2+ABS(B[j,i]);
    IF Tmp2>Tmp THEN Tmp:=Tmp2;
    END;
  Norme_Colonne:=Tmp
  END;

  PROCEDURE Init(N:BYTE;A:Matrice;VAR B:Matrice;VAR OK:BOOLEAN);
  VAR i,j:BYTE;
      T:REAL;
  BEGIN
  FILLCHAR(Id,SizeOf(Id),0);
  T:=Norme_Ligne(N,A)*Norme_Colonne(N,A);
  IF ABS(T)>1E-10 THEN BEGIN
                       OK:=TRUE;
                       FOR i:=1 TO N DO
                         BEGIN
                         Id[i,i]:=1;
                         B[i,i]:=A[i,i]/T;
                         FOR j:=i+1 TO N DO
                           BEGIN
                           B[i,j]:=A[j,i]/T;
                           B[j,i]:=A[i,j]/T;
                           END;
                         END;
                       END
  ELSE OK:=FALSE;
  END;

  PROCEDURE Calcule_E;
  VAR C:Matrice;
  BEGIN
  Produit(N,B,A,C);
  Difference(N,Id,C,E);
  END;

  PROCEDURE Calcule_B;
  VAR C,D:Matrice;
  BEGIN
  Somme(N,Id,E,D);
  Produit(N,D,B,C);
  B:=C
  END;

BEGIN
Init(N,A,B,OK);
IF NOT OK THEN WRITELN('Matrice nulle!')
ELSE BEGIN
     Calcule_E;
     Old:=Norme_Ligne(N,E);
     New:=Old;
     Cpt:=0;
     WHILE NOT ((New<1E-6) OR ((New>Old) AND (Cpt>5))) DO
       BEGIN
       Old:=New;
       Calcule_B;
       Calcule_E;
       New:=Norme_Ligne(N,E);
       INC(Cpt);
       END;
    OK:=TRUE;
    IF New<1E-6 THEN X:=B
    ELSE OK:=FALSE
    END;
END;

BEGIN
Intro(N,MatA);
Affiche('A=',N,MatA);
Attente;
Inverse(N,MatA,MatB,OK);
IF OK THEN BEGIN
           Affiche('Inverse de A=',N,MatB);
           Attente;
           Produit(N,MatA,MatB,MatC);
           Affiche('Produit des matrices inverses=',N,MatC);
           Attente;
           END
ELSE BEGIN
     CLRSCR;
     WRITELN('Pas d''inverse!');
     END;
END.


Cliquez ici pour charger une archive contenant les implémentations de l'algo en Turbo Pascal et en Delphi.

Page d'accueil.