Orthonormalisation de gram-schmidt

Description

Programmation du procédé d'orthonormalisation de Gram-Schmidt avec des vecteurs de R^n, et avec un produit scalaire quelconque.

Source / Exemple :


#ifndef _UTIL_H_
	#include "util.h"
#endif // _UTIL_H_

#ifndef _MATH_H_
	#include "math.h"
#endif // _MATH_H_

#ifndef _LIST_H_
	#include "list.h"
#endif // _LIST_H_

#ifndef _WINUTIL_H_
	#include "winutil.h"
#endif // _WINUTIL_H_

#include "conio.h"

//-------------------------------------------------
// OUTILS
//-------------------------------------------------
void printFloat(double x)
{
printf("%lf",x);
} // printFloat()

//-------------------------------------------------
void printVect(P_VECTND u)
{
int i;
printf("(");
printFloat(u->coord[0]);
for(i=1;i<u->dim;i++)
  {
  printf(" ,");
  printFloat(u->coord[i]);
  }
printf(" )\n");
} // printVect()

//-------------------------------------------------
// GRAM-SCHMIDT
//-------------------------------------------------
typedef double (*DOT_PRODUCT)(P_VECTND,P_VECTND);

//-------------------------------------------------
double dotProduct(DOT_PRODUCT phi,P_VECTND u,P_VECTND v)
{
return phi(u,v);
} // dotProduct()

//-------------------------------------------------
double norm(DOT_PRODUCT phi,P_VECTND u)
{
return sqrt(phi(u,u));
} // norm()

//-------------------------------------------------
PP_VECTND GramSchmidt(PP_VECTND vect,int nb,DOT_PRODUCT phi)
{
PP_VECTND res;
int       dim,i;
double    n;

AssertPointer(vect);
AssertPointer(phi);
Assert(nb >= 1);

dim = vect[0]->dim;
res = Malloc(P_VECTND,nb);
for(i=0;i<nb;i++)
  {
  res[i] = CreateVectND(dim);
  }

// le premier vecteur est juste normalise
n = norm(phi,vect[0]);
Assert(!CmpDouble0(n));
MulByScalarVectND(res[0],vect[0],1./n);
for(i=1;i<nb;i++)
  {
  int j;

  Assert(vect[i]->dim == dim);
  CopyVectND(res[i],vect[i]);

  // on hote les composantes suivants les vecteurs deja "orthonormalises"
  for(j=0;j<i;j++)
    {
    double k;
    k = dotProduct(phi,res[i],res[j]); // la composante est le produit scalaire ...
    LinCombinVectND(res[i],res[i],res[j],1.,-k);
    }
  
  // on normalise
  n = norm(phi,res[i]);
  Assert(!CmpDouble0(n));
  MulByScalarVectND(res[i],res[i],1./n);
  }

return res;
} //GramSchmidt()

//-------------------------------------------------
// MAIN
//-------------------------------------------------
double myDotProduct(P_VECTND u,P_VECTND v)
{
double x1,y1,z1;
double x2,y2,z2;

x1 = u->coord[0];
y1 = u->coord[1];
z1 = u->coord[2];

x2 = v->coord[0];
y2 = v->coord[1];
z2 = v->coord[2];

return x1*x2+y1*y2+z1*z2+0.5*(x1*y2+x2*y1+x1*z2+x2*z1+y1*z2+y2*z1);
} // myDotProduct()

//-------------------------------------------------
int main(int argc,char **argv)
{
InitializationLibUtil();
  {
  P_VECTND  vect[3];
  PP_VECTND res;

  vect[0] = CreateVectND(3);vect[0]->coord[0] = 1;vect[0]->coord[1] = -1;vect[0]->coord[2] = 0;
  vect[1] = CreateVectND(3);vect[1]->coord[0] = 1;vect[1]->coord[1] =  1;vect[1]->coord[2] = 0;
  vect[2] = CreateVectND(3);vect[2]->coord[0] = 1;vect[2]->coord[1] = -1;vect[2]->coord[2] = 1;

  res = GramSchmidt(vect,3,myDotProduct);

  printVect(vect[0]);
  printVect(vect[1]);
  printVect(vect[2]);
  printf("\nPar le procede d'orthonormalisation de Gram-Schmidt\n");
  printf("la famille de vecteurs ci-dessus devient :\n\n");
  printVect(res[0]);
  printVect(res[1]);
  printVect(res[2]);

  DeleteVectND(vect[0]);
  DeleteVectND(vect[1]);
  DeleteVectND(vect[2]);

  DeleteVectND(res[0]);
  DeleteVectND(res[1]);
  DeleteVectND(res[2]);

  Free(res);
  }
CloseLibUtil();
CheckingClosingLibUtil();

printf("\nfin du programme ...\n");
getch();
return 0 ;
} // main()

Codes Sources

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.