Polynome caracteristique

Soyez le premier à donner votre avis sur cette source.

Vue 17 189 fois - Téléchargée 575 fois

Description

- algorithme pour calculer le polynome caracteristique d'une matrice (carree) :
Pc(M)(x) = det(M - x.I)

- l'algorithme consiste a faire un FFT a partir de certaines valeurs du polynome
que l'on sait calculer avec le determinant d'un matrice (algo. en n^3)

- trois fonctions :
  • det -> calcule le determinant d'un matrice (complexe)
  • polyFFT -> calcule la FFT d'un tableau de nombres complexes
  • Pc -> retourne le polynome caracteristique


- voir captures d'ecran pour des exemples de matrices

Source / Exemple :


//-------------------------------------------------
// DET
//-------------------------------------------------
// calcul d'un determinant avec le pivot de Gauss
// la complexite est en n^3
COMPLEX *Det(COMPLEX **_mat,int n,COMPLEX *res)
{
COMPLEX **mat;
COMPLEX   det;
int       x,y;
int       i,j;

mat = CreateMat(n);
for(y=0;y<n;y++)
  {
  for(x=0;x<n;x++)
    {
    mat[y][x].re = _mat[y][x].re;
    mat[y][x].im = _mat[y][x].im;
    }
  }

det.re = 1.;
det.im = 0.;

for(j=0;j<n-1;j++)
  {

  int       rankMax,rank;
  COMPLEX   coeffMax;

  // ( etape 1 )
  rankMax = j;
  for(rank=j+1;rank<n;rank++)
    {
    if(AbsComplex(&mat[rankMax][j]) < AbsComplex(&mat[rank][j]))
      {
      rankMax = rank;
      }
    }

  coeffMax.re = mat[rankMax][j].re;
  coeffMax.im = mat[rankMax][j].im;
  if(AbsComplex(&coeffMax) <= MATH_EPSILON)
    {
    det.re = 0.;
    det.im = 0.;
    goto label_end;
    }
  // ( etape 2 )
  if(rankMax != j)
    {
    for(i=j;i<n;i++)
      {
      COMPLEX tmp;
      tmp             = mat[j][i];
      mat[j][i]       = mat[rankMax][i];
      mat[rankMax][i] = tmp;
      }
    det.re *= -1.;
    det.im *= -1.;
    }

  MulComplex(&det,&det,&coeffMax);
  // ( etape 3 )
  for(rank=j+1;rank<n;rank++)
    {
    COMPLEX coeff;
    DivComplex(&coeff,&mat[rank][j],&coeffMax);
    for(i=j;i<n;i++)
      {
      COMPLEX mul;
      MulComplex(&mul,&coeff,&mat[j][i]);
      SubComplex(&mat[rank][i],&mat[rank][i],&mul);
      }
    }
    
  }

MulComplex(&det,&det,&mat[n-1][n-1]);

label_end:
DeleteMat(mat,n);

  • res = det;
return res; } // Det() //------------------------------------------------- // POLY_FFT //------------------------------------------------- // on decompose le polynome P = A0 + B0 avec les coefficients // de degres pairs (en A0) et impairs (en B0) // alors on a : // P(w) = A0(w) + B0(w) // = A(w^2) + w.B(w^2) // or w^2 est une racine nieme, donc on se ramene a un cas // de taille deux fois moins grand void polyFFT( P_COMPLEX res, P_COMPLEX arrayCoeff, int stepi, int n, P_COMPLEX w ) { if(1 == n) { // P(z)=constante res->re = arrayCoeff->re; res->im = arrayCoeff->im; } else { int k; COMPLEX w2,wPowk; // strategie diviser pour regner // calcul du carre de <w> MulComplex(&w2,w,w); // appel recursif polyFFT(res ,arrayCoeff ,2*stepi,n/2,&w2); // coefficients pairs polyFFT(res + n/2,arrayCoeff + stepi ,2*stepi,n/2,&w2); // coefficients impairs // on calcule la FFT a partir de ce que l'on vient de calculer wPowk.re = 1.; wPowk.im = 0.; for(k=0;k<n/2;k++) { COMPLEX A,B; A = res[k]; B = res[k + n/2]; MulComplex(&B,&B,&wPowk); AddComplex(res + k ,&A,&B); SubComplex(res + k + n/2,&A,&B); MulComplex(&wPowk,&wPowk,w); } } } // polyFFT() //------------------------------------------------- // PC //------------------------------------------------- // // calcul du polynome caracteristique // Pc(M)(x) = det(M - x.I) // P_POLYNOMIAL Pc(COMPLEX **_mat,int n) { int len; // aligne sur une puissance de 2 P_POLYNOMIAL res; // le polynome caracteristique COMPLEX *arrayDet; // valeurs des determinants COMPLEX **mat; int k,x,y; COMPLEX conjw; len = alignedPower2(1 + n); res = CreatePolynomial(len - 1); arrayDet = Malloc(COMPLEX,len); mat = CreateMat(n); // on copie la matrice for(y=0;y<n;y++) { for(x=0;x<n;x++) { mat[y][x].re = _mat[y][x].re; mat[y][x].im = _mat[y][x].im; } } // on evalue le polynome aux racines <len>-eme de l'unite for(k=0;k<len;k++) { COMPLEX w; int i; GetNthRootUnitComplex(&w,k,len); for(i=0;i<n;i++) { SubComplex(&mat[i][i],&_mat[i][i],&w); } // calcul P(w) dans <arrayDet + k> Det(mat,n,arrayDet + k); } // on effectue une FFT pour retrouver les coefficients du polynome // a partir de ses valeurs en certains points GetNthRootUnitComplex(&conjw,-1,len); polyFFT(res->coeff,arrayDet,1,len,&conjw); CalcDeg(res); for(k=0;k<=res->deg;k++) { res->coeff[k].re /= len; res->coeff[k].im /= len; } Free(arrayDet); DeleteMat(mat,n); return res; } // Pc()

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.