TPs Calcul Numérique

VI.1 - Méthode De Krylov

Objectif : calcul des coefficients du polynôme caractéristique

Cette méthode repose sur la résolution du système linéaire BX=C Où la matrice B est construite comme suit :

On note les colonnes de la matrice B par bj ; j : 1, n Soit X0 vecteur donné.

bn=X0

bn-1=AX0 bn-2=A²X0

.
.
.

b1 =Aⁿ+¹X0

et le vecteur C=AⁿX0

et par l'utilisation de la méthode de Gauss on obtient n coefficients + le premier (+1) n.

Le programme

//TP10 Calcul Numérique

//Calcul des coefficients du polynome caracteristique //Méthode de Krylov

#include<iostream.h> #include<iomanip.h> #include<conio.h> #include<math.h> #include<alloc.h>

//////////////////Structure pour le programme de Gauss///////////////////////

struct PosMax //enregistrement pour le maximum

{int IL;int JC;};

//enregistrement pour la solution

struct solution {double x;int indice;};

//////////////////////////////Fin////////////////////////////////////////////////

double **DefMat(int n,int m)//allocation dynamique d'une matrice {double * *M=(double* *)malloc(n*sizeof(double*));

for(int i=0;i<n;i++)

M[i]=(double*)malloc(m*sizeof(double));

return M;

}

void remplir(double **A,int n,int m)//remplissage d'une matrice

{int x=0,y=wherey(); for(int i=0;i<n;i++)

{for(int j=0;j<m;j++)

{gotoxy(x+=8,y);

cin>>A[i] [j]; }

x=0;

y++;

}

}

void print(double **A,int n,int m)//affichage d'une matrice

{int x=0,y=wherey(); for(int i=0;i<n;i++)

{for(int j=0;j<m;j++)

{gotoxy(x+=8,y);

cout<<setw(8)<<setfill(' ')<<setprecision(2)<<A[i] [j]; }

x=0;

y++;

}

}

//produit d'une matrice A avec vecteur V

double *ProduitAV(double **A,double *V,int n)

{double *tmp=(double*)malloc(n*sizeof(double));

for(int i=0;i<n;i++)
{double som=0;

for(int l=0;l<n;l++) som+=A[i] [l] *V[l];

tmp[i]=som;

}

return tmp;

}

//calcul de A puissance P fois le vecteur V

double *puissance(double **A,double *V,int n,int p)

{

for(int i=0;i<p;i++) V=ProduitAV(A,V,n);

return V;

}

//-1 a la puissance de n int MoinsUn(int n) {int tmp=-1;

for(int i=1 ;i<n;i++) tmp*=-1;

return tmp;

}

//////////////////////////Methode de Gauss//////////////////////////////////// //Chercher la posision du max

PosMax SearchMax(double **mat,int n,int m,int k)

{PosMax max;

max.IL=max.JC=k;

double tmp=fabs(mat[k] [k]);

for(int i=k;i<n;i++) for(int j=k;j<m;j++)

if(tmp<fabs(mat[i] [j]))

{max.IL=i;max.JC=j;

tmp=fabs(mat[i] [j]);

}

return max;

}

//fonction de la permutation

void permut(double **mat,int n,int m,int k,solution *s)

{PosMax grand=SearchMax(mat,n,n,k);

if(grand.IL != k)

{double inter;

for(int j=k;j<m;j++)

{inter=mat[k] [j] ;mat[k] [j]=mat[grand.IL] [j] ;mat[grand.IL] [j]=inter; }

}

if(grand.JC != k)

{double inter;

for(int i=k;i<n;i++)

{inter=mat[i] [k] ;mat[i] [k]=mat[i] [grand.JC] ;mat[i] [grand.JC]=inter; }

int ind=s[k] .indice;s[k] .indice=s[grand.JC] .indice;s[grand .JC] .indice=ind; }

}

//fonction de la triangularisation

int triangul(double **mat,int n,int m,solution *s,double eps=0.00000000001) {for(int k=0;k<n-1 ;k++)

{permut(mat,n,m,k,s);

if(fabs(mat[k][k]) > eps)

for(int i=k+1 ;i<n;i++)

{double pivot=(double)mat[i] [k]/(double)mat[k] [k];

for(int j=k+1 ;j<m;j++) mat[i] [j]-=pivot*mat[k] [j];

mat[i] [k]=0;

}

else {cout<<"Matrice Singuliere ! ! ! ";return 0; }

}

if(fabs(mat[n- 1] [n-1]) <= eps) {cout<<"Matrice Singuliere ! ! ! ";return 0; } return 1;

}

//fonction de resolution

void solve(double **mat,int n,int m,solution *s)

{int t=n-1,l=m-1;

s [t] .x=(double)mat[t] [l]/(double)mat[t] [t];

for(int i=t-1 ;i>=0;i--)

{double som=0.0;

for(int j=i+1 ;j<n;j++) som+=mat[i] [j] *s [j] .x;

s[i] .x=(double)(mat[i] [l] - som)/(double)mat[i] [i];

}

}

//mis on ordre les solution

void ordre(solution *s,int n)

{

for(int i=0;i<n;i++)

{int ind=s[i] .indice;

for(int j=i+1 ;j<n;j++) if(ind>s[j] .indice) ind=s[j] .indice;

if(ind != i)

{solution tmp=s[i] ;s[i]=s[ind] ;s[ind]=tmp; }

}

}

//test des solutions

int test(double **A,int n,int m,solution *s,double eps=0.3)

{

for(int i=0;i<n;i++)

{double som=0;

for(int j=0;j<n;j++) som+=A[i] [j] *s[j] .x;

double min=A[i] [m-1] -eps,max=A[i] [m-1 ]+eps; if((som<min)||(som>max)) return 0;

}

return 1;

}

/////////////////////////////////////Fin////////////////////////////////////// void main()

{clrscr();

cout<<"\t\tMethode De Krylov\n";

cout<<"Donner Le Dimension De La Matrice :";

int n;cin>>n;//dimension de la matrice

double * *C=DefMat(n,n);

cout<<"Donner La Matrice \n";

remplir(C,n,n);//remplissage de la matrice

double *V=(double*)malloc(n*sizeof(double)); cout<<"Donner Le Vecteur V \n";

for(int j=0;j<n;j++) cin>>V[j] ;//la saisie du vecteur initial double **B=DefMat(n,n);//allocation dynamique de la mtrice B //construction de la matrice B

for(int i=0;i<n;i++)

{double *coef=puissance(C,V,n,i);

for(j=0;j<n;j++) B [j] [n-i-1 ]=coef[j];

}

//affichage de la matrice B

cout<<"La matrice B\n";

print(B,n,n);

cout<<endl;

//construction du vecteur Xn

double *Xn=puissance(C,V,n,n);

cout<<"Le vecteur Xn"<<endl;

//affichage du vecteur Xn

for(j=0;j<n;j++) cout<<Xn[j]<<endl;

////////////////////////////////Methode de Gauss///////////////////////////////// int dim;

dim=n;

solution *X=(solution*)malloc(dim*sizeof(solution)); for(i=0;i<dim;i++) X[i] .indice=i;//initialisation des indices

//Declaration de la matrice augmente double * *A=DefMat(dim,dim+1); //initialisation de A

for(i=0;i<n;i++)

for(j=0;j<n;j++)

A[i][j]=B[i][j];

for(i=0;i<n;i++) A[i] [n]=Xn[i];

int ok=triangul(A,dim,dim+1 ,X);//triangularisation if(ok)

{ cout<<"Matrice Triangulée\n";

print(A,dim,dim+1 );//A triangule

cout<<endl;

solve(A,dim,dim+1 ,X);//resolution

ordre(X,dim);//ordonancement des solutions ok=test(A,dim,dim+1 ,X);//teste des solutions if(ok)

{ cout<<"Les Coefficient Du Polynome Caracteristique Sont :\n"; cout<<endl;

for(i=0;i<n;i++) X[i] .x*=(-1 );//multiplication des solutions par -1

int t=MoinsUn(dim);//multiplication des solutions par -1 puissance n for(i=0;i<n;i++) X[i] .x*=t;

dim%2?cout<<"a0 = -1 ":cout<<"a0 = 1 ";cout<<endl;//affichage du 1er coef for(i=0;i<dim;i++) cout<<"a"<<(X[i] .indice)+1 <<" = "<<(X[i] .x)<<endl;

}

else cout<<"La solution a une erreur plus grande"; }

free(X);free(A);//liberation memoire

////////////////////////////////Fin//////////////////////////////////////////// free(C);free(V);free(B);free(Xn);//liberation memoire getch();

}