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=A2X0
.
.
.
b1 =An+1X0
et le vecteur C=AnX0
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();
}
| 
