Quadriques et reduction de gauss

5/5 (6 avis)

Vue 13 011 fois - Téléchargée 593 fois

Description

J'ai ete tres supris de ne pas voir sur ce site un code sur la reduction des formes quadratiques
reelles par l'algorithme de reduction de Gauss ...

Les quadratiques sont l'equivalent des coniques en dimensions 3, et il existe une classification
des quadriques (ce sont des nappes, des surfaces) avec des noms affreux !

Pour classifier les quadriques on a besoin de la signature de leurs formes quadratiques associees
(il y en a deux notee M et M0). La signature d'un f.q.r. (forme quadratique reelle) est un couple (p,q)
de deux entiers : si la f.q.r. est mise sous forme exclusivement de carre (par diagonalisation
de la matrice symetrique reelle par exemple, ou par reduction de Gauss (cf. matrice congruente
a une matrice diagonale)) alors le premier entier est le nombre de coefficients strictements
positifs et le second pour les negatifs. La signature d'un f.q.r. ne depend pas de la base
d'expression, i.e. elle est invariante pour deux matrices (symetriques) congruentes.
Je rappelle que la congruence entre A et B s'exprime par A=tP.B.P avec tP la transposee de P, et P
une matrice inversible, a ne surtout pas confondre avec la similitude A=P^-1.B.P, alors
les endomorphismes sont identiques.

D'un point de vu mathematique, la reduction de Gauss permet de classifier les classes de
congruences, ceci est en rapport direct avec la loi d'inertie de Sylvester dans le cas ou le
corps de base est reel. Une autre demonstration de l'existence de la reduction (mais ne donne
pas d'algorithme pour la trouver, a part si l'on sait diagonaliser une matrice ... pas
si facile que cela) est le theoreme spectral : toute matrice symetrique reelle est diagonalisable
dans une base orthonormee (tout est dans le orthonormee ... P^-1=tP ...)

La reduction de Gauss est une methode d'elimination, on cherche a mettre notre f.q.r. sous la
forme d'une somme d'une term carre plus une autre f.q.r. dependant d'un parametre de moins,
en fait il y a certains cas ou l'on mets sous la forme de 3 termes : deux carres et une f.q.r.
avec deux parametres en moins.

Source / Exemple :


//-------------------------------------------------
// GAUSS
//-------------------------------------------------
void gauss(double **M,int N)
{
int n,x,y,k;

n = N-1;
while(n > 0)
  {

  // si on n'a pas un terme carre ...
  if(CmpDouble0(M[n][n]))
    {
    int s; // indice du "swap"

    // ... on cherche un autre terme carre ...
    s = -1;
    for(k=0;k<n;k++)
      {
      if(!CmpDouble0(M[k][k]))
        {
        s = k;
        break;
        }
      }
    
    if(-1 == s)
      {
      // ... si on n'en trouver decidement pas, on regarde si on a des produits
      // croise entre la derniere variable et l'avant derniere, puis on elimine
      // les deux dernieres lignes et colonne :
      if(!CmpDouble0(M[n][n-1]))
        {
        // 
        // ( *   *   *   0   0    )
        // ( *   *   *   0   0    )
        // ( *   *   *   0   0    )
        // ( 0   0   0  1/2  0    )
        // ( 0   0   0   0  -1/2  )
        // 
        for(y=0;y<n-1;y++)
          {
          for(x=0;x<n-1;x++)
            {
            M[y][x] -= (M[n-1][x]*M[n][y]+M[n-1][y]*M[n][x])/M[n][n-1];
            }
          }
        M[n-1][n-1] = + 0.5;
        M[n-1][n  ] =   0.0;
        M[n  ][n-1] =   0.0;
        M[n  ][n  ] = - 0.5;

        // on diminue de 2 le rang
        n -= 2;
        }
      else
        {
        // sinon
        // 
        // ( 0   *   *   *   0  )
        // ( *   0   *   *   0  )
        // ( *   *   0   *   0  )
        // ( *   *   *   0   0  )
        // ( 0   0   0   0   0  )
        // 

        // on diminue de 1 le rang
        n -= 1;
        }
      }
    else
      {
      // ... si on en a trouve un, on l'echange et on recommence la boucle au meme rang
      // mais cette fois-ci avec un terme carre
      for(k=0;k<=n;k++) // ligne : L[n]  <->  L[s]
        {
        SwapDouble(&M[n][k],&M[s][k]);
        }
      for(k=0;k<=n;k++) // colonne : C[n]  <->  C[s]
        {
        SwapDouble(&M[k][n],&M[k][s]);
        }
      }
    }
  else
    {
    // ... on a donc un terme carre, on elimine la derniere ligne et colonne :
    // 
    // ( * * * * 0 )
    // ( * * * * 0 )
    // ( * * * * 0 )
    // ( * * * * 0 )
    // ( 0 0 0 0 ? )
    // 
    for(y=0;y<n;y++)
      {
      for(x=0;x<n;x++)
        {
        M[y][x] -= M[n][x] * M[y][n] / M[n][n];
        }
      }

    for(k=0;k<n;k++)
      {
      M[k][n] = M[n][k] = 0.;
      }

    // on passe au rang inferieur
    n --;
    }
  }
} // gauss()

//-------------------------------------------------
void sign(int *p,int *q,double **mat,int n)
{
double **diag;
int      k;

(*p) = 0;
(*q) = 0;

diag = CopyMat(mat,n);
gauss(diag,n);

for(k=0;k<n;k++)
  {
  if(diag[k][k] > MATH_EPSILON)
    {
    (*p) ++;
    }
  else if(diag[k][k] < -MATH_EPSILON)
    {
    (*q) ++;
    }
  }

DeleteMat(diag,n);
} // sign()

//-------------------------------------------------
void info(int pA,int qA,int pB,int qB)
{
if(4==pB && 0==qB && 3==pA && 0==qA)
  {
  printf("ensemble vide (quadrique propre)\n");
  }
else if(3==pB && 1==qB && 3==pA && 0==qA)
  {
  printf("ellipsoide (quadrique propre)\n");
  }
else if(3==pB && 1==qB && 2==pA && 1==qA)
  {
  printf("hyperboloide à deux nappes (quadrique propre)\n");
  }
else if(3==pB && 1==qB && 2==pA && 0==qA)
  {
  printf("paraboloide elliptique (quadrique propre)\n");
  }
else if(2==pB && 2==qB && 2==pA && 1==qA)
  {
  printf("hyperboloide à une nappe (quadrique propre)\n");
  }
else if(2==pB && 2==qB && 1==pA && 1==qA)
  {
  printf("paraboloide hyperbolique (quadrique propre)\n");
  }
else if(3==pB && 0==qB && 3==pA && 0==qA)
  {
  printf("point\n");
  }
else if(3==pB && 0==qB && 2==pA && 0==qA)
  {
  printf("ensemble vide\n");
  }
else if(2==pB && 1==qB && 2==pA && 1==qA)
  {
  printf("cone\n");
  }
else if(2==pB && 1==qB && 2==pA && 0==qA)
  {
  printf("cylindre elliptique\n");
  }
else if(2==pB && 1==qB && 1==pA && 1==qA)
  {
  printf("cylindre hyperbolique\n");
  }
else if(2==pB && 1==qB && 1==pA && 0==qA)
  {
  printf("cylindre parabolique\n");
  }
else if(2==pB && 0==qB && 2==pA && 0==qA)
  {
  printf("droite\n");
  }
else if(2==pB && 0==qB && 1==pA && 0==qA)
  {
  printf("ensemble vide\n");
  }
else if(1==pB && 1==qB && 1==pA && 1==qA)
  {
  printf("deux plans sécants\n");
  }
else if(1==pB && 1==qB && 1==pA && 0==qA)
  {
  printf("deux plans parallèles\n");
  }
else if(1==pB && 0==qB && 1==pA && 0==qA)
  {
  printf("plan double\n");
  }
else
  {
  printf("impossible\n");
  }
} // info()

Codes Sources

A voir également

Ajouter un commentaire Commentaires
Saros Messages postés 921 Date d'inscription vendredi 20 décembre 2002 Statut Membre Dernière intervention 23 septembre 2010
30 oct. 2007 à 20:03
A la main, ça se fait plutôt bien.. j'imagine que les applications doivent être assez spécifiques même si ça concerne de nombreux domaines (détermination d'extrema, équations différentielles ?)
En tout cas il t'a fallu du courage pour distinguer tous les cas :)
dletozeun Messages postés 546 Date d'inscription vendredi 13 février 2004 Statut Membre Dernière intervention 9 janvier 2008 1
5 août 2007 à 16:56
ok merci! ^^

Oui je suis tout a fait d'accord avec toi, ton code est vraiment très clair ( je voulais le préciser masi j'ai oublié ^^).
En tant qu'ancien taupin je confirme! Rien que pour calculer un polynome caractéristique correctement, on peut trouver 10 résultats différents! Et les quadriques, j'en ai plus aucun souvenir...^^
cs_JCDjcd Messages postés 1138 Date d'inscription mardi 10 juin 2003 Statut Membre Dernière intervention 25 janvier 2009 4
5 août 2007 à 16:39
C'est absolument pas brutal, et c'est tres correct :
en fait je me disais qu'il y avait beaucoup de code
sur l'inversion de matrices par pivot de Gauss, alors
qu'il existe d'autre algorithme comme celui-ci qui
permet de s'entrainer sur les indices, donc c'etait
pour montrer que l'on peut aisement diversifier...
et sinon je suis doit a fait d'accord que les f.q.r.
c'est complique pour un debutant en mathematique ...
j'y avais pas trop pense ... ;-)

Par contre je voulais aussi faire remarquer que mon
implementation n'est pas obscure, les indices sont
claires, par contre je suis d'accord que les disjonctions
de cas ne sont pas "lisibles" a premiere vue, bien que
l'indentation devrait aider ...

Ce programme peut aussi servir aux taupins qui ont
comme exercice d'application de leur cours sur les
quadriques de verifier la nature de leurs quadriques, car
quand on le fait a la main on se trompe toujours !
(je suis sur que l'on pourrait avoir plein de temoingnages)

Bien amicalement JCDjcd
dletozeun Messages postés 546 Date d'inscription vendredi 13 février 2004 Statut Membre Dernière intervention 9 janvier 2008 1
5 août 2007 à 12:33
Désolé je ne voulais pas te vexer... et ne me fais pas dire ce que j'ai pas dis s'il te plait.
"plus utile" ne veux pas dire utile...c'est juste que meme si ce code rentre parfaitement dans sa catégorie, j'ai vu des codes en faisant partie qui ont une réelle application dans le domaine de l'informatique et notamment de nombreux codes de ton cru. Tout ca pour te dire que je suis parfaitement d'accord avec le fait que les mathématiques sont énormément utiles en informatique, mais que ce code particulier a une application un peu trop spécifique à mon gout sans réelle utilité puisque je suis a peu pres certain que la plupart des logiciels de calcul formel savent calculer la signature d'une f.q.r....

Pour ce qui de l'"algorithme d'elimination tres bien pour les debutants qui veulent s'entrainer a programmer des petits algorithmes marrants, c'est tres instructif pour la manipulation des indices"

Je suis assez septique, encore faut t'il des maitriser le sujet des f.q.r et tout ce qui semble aller avec: diagonalisation, reductiuon de Gauss, forme exclusivement carrée...c'est plutot rebutant.
Je pense qu'une fois avoir compris tout ca on maitrise tres bien les manipulations d'indices dans une matrice, il suffirait simplement a un debutant de coder un algorithme de multiplication de matrices pour s'en faire une bonne idée...

Cordialement, en espérant ne pas paraitre brutal et incorrect. ^^
cs_JCDjcd Messages postés 1138 Date d'inscription mardi 10 juin 2003 Statut Membre Dernière intervention 25 janvier 2009 4
5 août 2007 à 00:12
Ha parceque tu crois qu'il y a beaucoup de codes utiles ??
Je pense (sns aucun pretention) que ce code est beaucoup
plus utile que la moitier des sources sur ce site pour
la simple raison que deja ce n'est pas un code sur un
truc deja fait 100 fois, puis c'est dans la rubrique
maths & algorithme donc c'est plus utile (... ca c'est
mon point de vue qui n'engage que moi ...), et surtout
parceque c'est un petit algorithme d'elimination tres bien
pour les debutants qui veulent s'entrainer a programmer
des petits algorithmes marrants, c'est tres instructif pour
la manipulation des indices, tout comme les methodes de
pivots d'ailleurs.

Cordialement JCDjcd (sans avoir voulu etre trop "pretentieux")

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.