Matrices réelles

Description

Bonjour,

Voici un ensemble élémentaire et traditionnel de routines de calcul matriciel réel.
Il comprend les opérations de base ainsi que le calcul du déterminant et de l'inversion, mais ne traite pas les valeurs propres et n'avantage pas les matrices particulières (diagonales, triangulaires, symétriques, creuses (ou éparses, en: sparse), etc).
Les vecteurs y sont considérés comme des matrices dont l'une des dimensions est 1 (vecteur 'ligne' ou vecteur 'colonne').

Les codes comparables (qui affichent la source) présentent parfois l'une ou l'autre des lacunes suivantes:
- calcul peu précis ou parfois même incorrect !
- généralisation de la règle de Sarrus qui est limitée à 3×3 !
- choix direct de l'élément diagonal comme pivot: blocage éventuel.
- choix du premier élément non nul comme pivot: imprécision.
- absence de tests de précision et de temps de calcul.

En 2011, j'ai déjà présenté un article similaire dans Javascript: Vecteurs et Matrices,
complété par Systèmes d'équations à n inconnues.

Essayez le programme interactif: Résolution de systèmes de N équations.

Remarque: normalement, on n'utilise pas l'inversion pour résoudre un système d'équations, car il existe des méthodes plus efficaces comme par exemple la décomposition en matrices triangulaires.
 
 

MatWV.h


Le fichier MatWV.h devrait vous donner une information suffisante pour pouvoir utiliser les routines proposées:
#include <stdint.h>
#include <string>

typedef unsigned int uint;
typedef double real; // or typedef float real;
using namespace std;

class MatWV {
public:
////// Constructeurs:
  MatWV(uint n);           // Matrice unité n×n
  MatWV(uint,uint);        // Matrice nulle nCol×nLin
  MatWV(uint,uint,real*);  // Matrice nCol×nLin, copie des éléments
  MatWV(uint,uint,real);   // Matrice aléatoire avec valeurs abs limitées
  ~MatWV();                                            // {conditions}
////// Méthodes qui modifient la matrice T (this):
  void Set(uint c,uint l,real e); // T[c,l] = e
  void Add(real r);        // T = T + r
  void Add(MatWV *M);      // T = T + M                    {same dims}
  void Mul(real r);        // T = T * r
  void Mul(MatWV *M);      // T = T * M                    {T.nC=M.nL}
  void Sub(MatWV *M);      // T = T - M                    {same dims}
  void Transp();           // T = transp(T)
  real Invert();           // T = invert(T); return det    {T.nC=T.nL}
  void Unit();             // T = T / Norm(T): Vector Unit {T.nC=1 or T.nL=1}
////// Méthodes qui ne modifient pas T (this):
  bool Dim(MatWV *M);      // dims(T) =? dims(M)
  uint nCol();             // nombre de colonnes
  uint nLin();             // nombre de lignes
  real Get(uint c,uint l); // return T[c,l]
  real MatWV::Det();       // Det(T)                       {T.nC=T.nL}
  real MatWV::Norm();      // Norm(T)                      {nCol==1 ou nLin==1}
  real MatWV::EcartI();    // Ecart abs max de T / matrice unité
  string Str(string,string); // MatWV to String(eol,sep)
////// Méthodes qui créent une matrice (new) sans modifier T (this):
  MatWV *MatAdd(real r);   // return (T + r)
  MatWV *MatAdd(MatWV *M); // return (T + M)               {same dims}
  MatWV *MatCopy();        // return Copy(T)
  MatWV *MatInvert();      // return inversion(T)          {T.nC=M.nL}
  MatWV *MatMul(real r);   // return (T × r)
  MatWV *MatMul(MatWV *M); // return (T × M)[M.nC × nL]    {T.nC=M.nL}
  MatWV *MatSub(MatWV *M); // return (T - M)               {same dims}
  MatWV *MatTransp();      // return transp(T)
////// Opérateurs:
  // Ils seront ajoutés ultérieurement
private:
  real *Ele;  // Eléments ligne par ligne
  uint nC,nL; // nombre de colonnes > 0, nombre de lignes > 0
};

real Determinant(uint n,real *m); // modifie m[n*n] = éléments ligne par ligne

 
 

Inversion


L'inversion de matrice proposée correspond à une version classique non récursive, avec les particularités suivantes:
- Une seule routine.
- On opère sur la matrice elle-même.
- Choix du plus grand pivot possible (Gauss-Jordan): améliore la précision.
- La boucle la plus intérieure est bien optimisée.
- On retourne le déterminant qui est calculé en même temps.

Pour estimer la précision de l'inversion, on multiplie la matrice par son inverse:
Plus ce produit est proche de la matrice unité, plus le calcul est précis !
La routine EcartI permet de quantifier la précision.
 
 

Main.cpp


Le programme Main.cpp comporte un ensemble de vérifications, ainsi que des tests de précision et du temps d'exécution.
Son output figure dans le fichier Resultats.txt.
Observez-bien ces résultats, et le cas échéant, adaptez Main.cpp à d'autres tests qui peuvent vous intéresser.

Avez-vous remarqué qu'avec de mêmes dimensions, l'inversion ne prends qu'environ 1.5 fois le temps d'une multiplication ?

Par exemple, pour des matrices 500×500 et 1000×1000, mon portable (I5-2.3GHz) utilise:
- 0.047 s, 0.53 s pour le calcul du déterminant.
- 0.302 s, 8.95 s pour une multiplication
- 0.577 s, 14.6 s pour l'inversion.

Autre constatation intéressante (dans Resultats.txt):
Pour les très grandes matrices examinées >= 800×800, malgré que le déterminant sorte du domaine de définition des réels 64 bits (double), l'inversion se fait correctement et avec une bonne précision ! Surprenant ...
 
 
Comme il est quasi impossible d'écrire un programme sans aucun bug, je vous serais reconnaissant de me signaler les erreurs constatées.
Merci d'avance.

Dites-moi si le développement de la fonction real Determinant(uint n,real *m); vous intéresse.
Il est effectivement basé sur 3 ou 4 règles "simples" du calcul matriciel.

Codes Sources

A voir également

Vous n'êtes pas encore membre ?

inscrivez-vous, c'est gratuit et ça prend moins d'une minute !

Les membres obtiennent plus de réponses que les utilisateurs anonymes.

Le fait d'être membre vous permet d'avoir un suivi détaillé de vos demandes et codes sources.

Le fait d'être membre vous permet d'avoir des options supplémentaires.