Quadriques et reduction de gauss

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

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.