|
République Algérienne Démocratique Et
Populaire
Centre Universitaire Larbi Ben M 'hidi
Oum El Bouaghi
Institut Des Sciences Exactes
Département D'Informatique
Module : Calcul Numérique
Année Universitaire 2007+2008
TPs Calcul Numérique
Ecrit En Langage C++
3ème Année Ingénieurs
Encadré Par L'enseignant : A.Boukrouma
Réalisé Par L'étudiant : Merazga
Salim
Préface
Ce support contient un ensemble de méthodes d'Analyse
Numérique programmées en langage C++.
Les programmes sont testés sur des exemples
théoriques et ont donné des résultats corrects.
Il serait important qu'ils soient testés avec des
valeurs expérimentales afin de mieux approcher la convergence et par
là optimiser le temps de programmation.
Il est nécéssaire de tenir compte que les
introductions aux méthodes sont restreintes et manquent beaucoup de
détailles. Pour plusieurs de détaille il faut revenir aux livres
d'Analyse Numérique.
Table De Matière
I - Résolution De
Systèmes D'équations
Linéaires page 04
I.1 - Méthode De Gauss Avec Pivot Total page
06
I.2 - Méthode De Gauss - Jordan page 10
(Utilisée pour calculer l'inverse d'une
matrice)
I.3 - Méthode De Gauss - Seidel page 14
II - Interpolation D'une Fonction
Par Un Polynôme page 17
II.1 - Méthode De Newton page 18
II.2 - Méthode De Lagrange page 20
III - Résolution D'une
Equation Non Linéaire De Type f(x) = 0
page 22
III.1 - Méthode De Newton - Rap hson page 22
III.2 - Accélérateur D'Aïtken page
24
IV - Analyse De Données
Expérimentales page 26
- Méthode Linéaire Des Moindres
Carrés
V - Calcul Approché Des
Eléments Propres (valeurs et vecteurs propres)
page 31
- Méthode De La Puissance
Itérée
VI - Calcul Approché Des
Coefficients D 'un Polynôme de Degré n
Page 35
VI. 1 - Méthode De Krylov page 36
VI.2 - Méthode De Sauriau - Faddev page 41
VII - Recherche De La Valeur Approchée De La Plus
Grande Racine
D 'un Polynôme de
Degré n page 44
VII. 1 - Méthode De Graeffe page 44
VII.2 - Méthode De Bernoulli page 47
I - Résolution De Systèmes
D'équations Linéaires
En mathématique et plus exactement en algèbre
linéaire un système d'équation linéaire est un
ensemble d'équation linéaire, par exemple
13xi + 5x2 -- x3 = 2
2
5 xi + 3 x2 + x3 = 0
xi + 4x2 -- 9x3 = 1
Le problème qui se pose ici est comment peut-on
trouver x1, x2 et x3 qui satisfassent les trois équations
simultanément ?
Pour arriver à ce but on a plusieurs
méthodes, parmi ces méthodes on a la méthode de Gauss avec
pivot total, l'élimination de Gauss - Jordan et la méthode
itérative de Gauss - Seidel.
Généralement un système linéaire
peut être écrit sous la forme suivante :
aiixi + ai2x2 +
·
·
· + ainxn =
bi
Tel que
x j ; j : 1, n sont les inconnus
aij ; i : 1, m ; j : 1 ; n sont les coefficients
du système
{
anxi + a22x2 +
·
·
· + a2nxn =
b2
ami xi+am2 x2+
·
·
· + amnxn = bm
Ou sous la forme matricielle AX = B où
aii ai2
·
·
· ain xi
bi
a~i a22
·
·
· a2n
A =
·
· .X=x
·2et
B = b2 ~
ami am2
·
·
· amn xn
bm
Remarque : toutes les méthodes sont appliquées
sur des matrices carrées (nombre de ligne = nombre de colonne) ,sinon il
faut la rendre carrée par la multiplication du système au
transposé de la matrice elle-même
Supposons A une matrice non-carrée de n lignes et m
colonnes (de type (n,m)) et le système linéaire AX = B ;alors on
va faire le produit matricielle comme suit : At A X = At
B
Après le produit on obtient le nouveau systeme
linéaire W X = V tel que : W = At A (matrice carrée de
type (m, m) (m lignes et m colonnes))
V = At B (nouveau verteur de m
composantes)
La matrice At est obtenue par l'iversion des
lignes de la matrice A en des colonnes.
I.1 - Méthode De Gauss Avec Pivot Total Objectif
: résolution d'un système linéaire AX=B
On considère le système linéaire AX = B
suivant :
(Ill au '
·
· aln) (x1) (bi
a
n
1
an2
·
·
·
annxn
. . . .
.
. . . .
. . . . .
b2bn
=
Le but de cette méthode est de rendre la matrice A
triangulaire (supérieure ou inferieure) c'est-à-dire :
Matrice triangulaire supérieure
(a11 a12 '
·
· aln) (x1) (bi
0 a22
·
·
· a2n x2
. . .
. .. :
. :
. = b2
00
·
·
·
annxnbn
Et dans ce cas on peut trouver toutes (es valeurs de Xi par
la méthode itérative suivante (substitution
arrière)
ann
On carcure tout abord xn = bn
Puis
|
xi =
|
[1ii
|
n
bi -- 1 auxj ; i: n -- 1,1
j=i+1
|
Matrice triangulaire inferieure
(an 0 ... 0 ) (xi) (bi
a21 a22 ... 0 x2
. . . . .
. . . . . .
. . . . . .
_ b2
an
1
an2
·
·
·
annxnbn
De manière similaire on peut trouver toutes les
valeurs de Xi par la méthode itérative suivante
(substitution avant)
a11
On carcure tout abord xi = b1
Puis
)6
b -- + ajjXj / ;i:2-- 1,1
,-~
1
*
~))
X) ~
Le programme
//TP1 Calcul Numérique //Méthode de Gauss
#include<stdio.h> #include<alloc.h>
#include<iostream.h> #include<iomanip.h> #include<math.h>
#include<conio.h>
struct PosMax //enregistrement pour le maximum de la matrice {int
IL;int JC;};
//enregistrement pour la solution
struct solution {double x;int indice;};
double **NewMat(int n,int m) //creation dynamique d'une matrice
{double * *mat=(double* *)malloc(n*sizeof(double*));
for(int i=0;i<n;i++)
mat[i]=(double*)malloc(m*sizeof(double));
return mat;
}
void charger(double **mat,int n,int m) //initialisation du
matrice
{ cout<<"Entrer la matrice augmentée\n";
for(int i=0;i<n;i++)
{int x=wherex();
for(int j=0;j<m;j++)
{cin>>mat[i] [j] ;gotoxy(x+=6,wherey()- 1); }
cout<<endl;
x=8;
}
}
void print(double **mat,int n,int m)//afficher une matrice
{for(int i=0;i<n;i++)
{for(int j=0;j<m;j++)
cout<<setfill('
')<<setw(4)<<setprecision(2)<<mat[i] [j]<<"\t";
printf("\n");
}
}
//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);//recherche du max
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.0001)
{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;
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;
}
void main()
{char rep;
do
{clrscr();
cout<<"\t\t Resolution Du Systeme Linéaire
AX=B\n\t\t\tMethode de Gauss\n"; cout<<"Entrer la dimension du matrice
n:";
int dim;
cin>>dim;
solution *X=(solution*)malloc(dim*sizeof(solution));
for(int i=0;i<dim;i++) X[i] .indice=i;//initialisation des
indices
cout<<endl;
//Declaration de la matrice augmentée double
**A=NewMat(dim,dim+1);
charger(A,dim,dim+1 );//initialisation de A
int ok = triangul(A,dim,dim+1 ,X);//triangularisation
cout<<endl;
if(ok)
{ cout<<"Matrice Triangulée\n";
print(A,dim,dim+1 );//A triangulée
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 Solutions Sont :\n";
cout<<endl;
for(i=0;i<dim;i++) cout<<"X"<<X[i]
.indice<<"="<<X[i] .x<<endl;
}
else cout<<"La solution a une erreur plus grande"; }
free(A);free(X);
cout<<"Voulez vous répéter o/n:";
cin>>rep;
free(A);free(X);
}while (rep == 'o');
}
I.2 - Méthode De Gauss - Jordan
(Utilisée pour calculer l'inverse d'une matrice)
Objectif : calcul de la matrice inverse
On considère le système linéaire AX = B
suivant :
(a11 a12
·
·
· aln) (x1) (b1
am a22
·
·
· a2n
. ... .
x2
. = b2 .
an1 an2
·
·
· ann
xn
bn
1 0 0
(01 ... 0) x1 b1( bi
b2
=
bn
00
·
·
·1xn
Le but de cette méthode est de rendre la matrice A
diagonale unitaire c'est-à- dire :
. .
· . . .
. . . . . .
. . . . .
On peut trouver toutes les valeurs de Xi par la
méthode itérative suivante
xi = bi; i: 1,n
L'un des grands défis de l'algèbre
linéaire est de trouver la matrice inverse A-1 .Avec cette
méthode on peut calculer la matrice inverse sans utiliser la
cof acte ur (a i j)
formule théorique (A_1)i] =det (A)
qui est impraticable, par la
résolution de n système au même temps A
A-1 = I (I matrice identité
(a11 a12
·
·
· aln 1 0 ... 0
am a22
·
·
· a2n A
_i= ~0 1
·
·
· 0
an1 an2
·
·
· ann 0 0
·
·
· 1
. ... .
·
·
·
·
·
·
·
·
. . . .
Le programme
//TP2 Calcul Numérique
//Méthode de Jordan Pour Calculer L'inverse d' une
matrice
#include<stdio.h> #include<alloc.h>
#include<iostream.h> #include<iomanip.h> #include<math.h>
#include<conio.h>
struct PosMax //enregistrement pour le maximum {int IL;int
JC;};
//enregistrement pour la matrice inverse struct inverse {double
*x;int indice;};
double **NewMat(int n,int m) //creation dynamique d'une matrice
{double * *mat=(double* *)malloc(n*sizeof(double*));
for(int i=0;i<n;i++)
mat[i]=(double*)malloc(m*sizeof(double));
return mat;
}
void charger(double **mat,int n,int m) //initialisation du
matrice { cout<<"Entrer la matrice \n";
for(int i=0;i<n;i++)
{int x=wherex();
for(int j=0;j<m;j++)
{cin>>mat[i] [j] ;gotoxy(x+=6,wherey()- 1); }
cout<<endl;
}
}
void print(double **mat,int n,int m)//afficher une matrice
{for(int i=0;i<n;i++)
{for(int j=0;j<m;j++)
cout<<setfill('
')<<setw(4)<<setprecision(2)<<mat[i] [j]<<"\t";
printf("\n");
}
}
//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 normalisation
void normalis(double **mat,int m,int k)
{double p=mat[k] [k];
for(int j=k;j<m;j++) mat[k] [j]=(double)mat[k]
[j]/(double)p;
}
//fonction de la permutation
void permut(double **mat,int n,int m,int k,inverse *inv)
{PosMax grand=SearchMax(mat,n,n,k);
if(grand.IL != k)//permutation de ligne
{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)//permutation de colonne
{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=inv[k] .indice;inv[k] .indice=inv[grand.JC]
.indice;inv[grand .JC].indice=ind; }
}
//fonction de reduction
int reduction(double **mat,int n,int m,inverse *inv,double
eps=0.001)
{for(int k=0;k<n;k++) {permut(mat,n,m,k,inv);
if(fabs(mat[k][k]) > eps)
{normalis(mat,m,k);//normalisation de la ligne k
for(int i=0;i<n;i++)
if( i != k)
{for(int j=k+1 ;j<m;j++) mat[i] [j] -=(double)mat[i] [k]
*mat[k] [j];
mat[i] [k]=0; }
}
else {cout<<"Matrice Singuliere ! ! ! ";return 0; }
}
return 1;
}
//fonction pour extraire la matrice invese
void GetInv(double **mat,int n,inverse *inv)
{for(int i=0;i<n;i++)
{inv[i] .x=(double*)malloc(n*sizeof(double));
for(int j=0;j<n;j++) inv[i] .x[j]=mat[i] [n+j];
}
}
//fonction d'affichge pour la matrice inverse
void PrintInv(inverse *inv,int n)
{for(int i=0;i<n;i++) {for(int j=0;j<n;j++)
cout<<setfill('
')<<setw(4)<<setprecision(2)<<inv[i] .x[j]<<"\t";
printf("\n");
}
}
//mise on ordre les ligne de la matrice inverse
void ordre(inverse *inv,int n)
{for(int i=0;i<n;i++)
{int ind=inv[i] .indice;
for(int j=i+1 ;j<n;j++) if(ind>inv[j] .indice) ind=inv[j]
.indice; if(ind != i)
{inverse tmp=inv[i] ;inv[i]=inv[ind] ;inv[ind]=tmp; }
}
}
//chargement automatique de la matrice identitée
void ChargerI(double **mat,int n,int m)
{for(int i=0;i<n;i++)
for(int j=n;j<m;j++)
if(i == j-n) mat[i][j]=1;
else mat[i][j]=0;
}
void main()
{char rep;
do
{clrscr();
cout<<"\t\t Calcul Du Matrice Inverse\n\t\t Methode de
Jordan - Gauss\n"; cout<<"Entrer la dimension du matrice n:";
int dim;
cin>>dim;
inverse *inv=(inverse*)malloc(dim*sizeof(inverse));
for(int i=0;i<dim;i++) inv[i] .indice=i;//initialisation des
indices
cout<<endl;
//Declaration de la matrice augmentée
double * *A=NewMat(dim,dim*2);
charger(A,dim,dim);//initialisation de A
ChargerI(A,dim,dim*2);//ajoute la matrice identitée
cout<<"\nMatrice Augmentée Non
Traitée"<<endl;
print(A,dim,dim*2);
int ok = reduction(A,dim,dim*2,inv);//reduction
cout<<endl; if(ok)
{ cout<<"Matrice Augmentée
Traitée"<<endl;
print(A,dim,dim*2);
GetInv(A,dim,inv);//Donner l'inverse dans inv
cout<<endl;
ordre(inv,dim);//mise on ordre inverse
cout<<"La Matrice Inverse"<<endl;
PrintInv(inv,dim);//afficher la matrice inverse
}
free(A);free(inv);//libération mmoire
cout<<"Voulez vous répéter o/n:";
cin>>rep;
}while (rep == 'o');
I.3 - Méthode De Gauss - Seidel
Objectif : Résolution d'un système
linéaire d'une manière itérative
On utilise cette méthode Pour les très grosses
matrices, Elle converge vers la solution progressivement et plus
rapidement.
Idée de base est de résoudre une
équation de type : F(x)=X ; avec une suite de type : Xk+1
=F(Xk) ;X0 donnée.
Si |Xk+1 - Xk | <î alors Xk
solution.
Pour faire ça on décompose la matrice A tel que
A = F - E ; alors
AX=BF(F--E)X=BIFX--EX=BIFX=EX+B
I X=F6'EX+F6'BIX=?X+K
Où
F=(
0
0 a22
~ ~ ~ ~
0 0
·
·
· ann
a11 0
~
~
0
1 S 0 0
a11
V
0 1 a22
S ~ 0 U
~ ~ ~ ~ U
0 0
·
·
· 1 ann
S T
alors F6' =
n
2n
H = ~0 --a2 ~ --a --a2 0 ~ --a
~ ~ ~ ~
--ani --an2 ~ 0
Donc à chaque itération on a :
t6 n
1
t 7W.8 = 7'.' 8 -- > atix, 7W8
X bt -- > at,x, Y ; j: 1, 2
attj=1 j=t+1
Remarque : pour assurer la convergence il faut que A soit
dominante. Le programme
//TP3 Calcul Numérique //Méthode Gauss-Seidel
#include<iostream.h> #include<iomanip.h>
#include<math.h> #include<conio.h> #include<alloc.h>
double **DefMat(int n,int m)//Allocation dynamique du matrice
{double * *tmp=(double* *)malloc(n*sizeof(double*)); for(int i=0;i<n;i++)
tmp[i]=(double*)malloc(m*sizeof(double));
return tmp;
}
void charger(double **M,int n,int m)//Chargement du matrice {int
x=0;
for(int i=0;i<n;i++)
{for(int j=0;j<m;j++)
{cin>>M[i] [j];
gotoxy(x+=8,wherey()- 1);
}
cout<<endl;
x=0;
}
}
int dominante(double **A,int n)//vérification du
domination {for(int i=0;i<n;i++)
{double som=0;
for(int j=0;j<n,j !=i;j++) som+=fabs(A[i] [j]);
if(som >= fabs(A[i][i]) ) return 0;
}
return 1;
}
double* Xkplus1(double **A,double *X0,int n)//calcul du X(k+1)
{double *tmp=(double*)malloc(n*sizeof(double));
for(int i=0;i<n;i++)
{double S 1=0,S2=0;
for(int j=0;j<=i-1 ;j++) S1+=A[i][j]*tmp[j]; for( j=i+1;j<n
;j++) S2+=A[i][j]*X0[j];
tmp[i]=(double)(A[i] [n] - S1 - S2)/(double)A[i] [i]; }
return tmp;
}
//test de convergence
short test(double *X0,double *eps,double *Xk,int n) {for(int
i=0;i<n;i++)
if( fabs(Xk[i] - X0[i]) > eps[i]) return 0;
return 1;
}
//resolution
double *Solve(double **A,double *X0,double *eps,int n,int kmax=6)
{for(int k=0;k<kmax;k++)
{double *Xk=Xkplus1 (A,X0,n);
short ok=test(X0,eps,Xk,n);
if(ok) return Xk;
X0=Xkplus 1 (A,X0,n);
}
return NULL;//en cas pas de solution en kiem itération
}
void main()
{char rep;
do
{clrscr();
cout<<" Méthode De Gauss - Seidel\n";
int dim;
cout<<"\nEntere La Dimension Du Matrice :";
cin>>dim;
//X0 initial
double *X0=(double*)malloc(dim*sizeof(double)); for(int
t=0;t<dim;t++) X0[t]=1;
//le vecteur epselon
double *eps=(double*)malloc(dim*sizeof(double)); for(
t=0;t<dim;t++) eps[t]=0.2;
//La matrice augmentée A
double * *A=DefMat(dim,dim+1);
charger(A,dim,dim+1 );//Chargement de A
cout<<endl;
int ok=dominante(A,dim);
if(ok)//test de la domination
{//Resolution
double * Sol=Solve(A,X0,eps,dim);
if(Sol) for(t=0;t<dim;t++) cout<<"X"<<t<<" =
"<<setprecision(2)<<Sol[t]<<endl; else cout<<"On a pas
une solution par cette méthode ! ! !";
free(Sol);//libération du mmoire
}
else cout<<"\nLa Matrice n'est pas dominante ! ! !\n";
//libération du mémoire
free(A);free(X0);free(eps);
cout<<"\nVoulez vous répéter o/n :
";cin>>rep;
}while (rep == 'o');
}
II - Interpolation D'une Fonction
Par Un Polynôme
On a un ensemble de n+1 points : x0 .... Xn et
pour chaque point on connaît la valeur d'une fonction f inconnue
(à-priori): f(x0) ... . f(xn)
La question qui se pose c'est : Quelle est la valeur de f sur
les points intermédiaire ?
Pour cela on doit supposer un modèle
mathématique de f qui est souvent un polynôme.
Pour l'interpolation polynomiale on a 2 méthodes
:
Méthode De Newton
Méthode De Lagrange
II.1 - Méthode De Newton
Objectif : Obtention d'un polynôme à
partir d'une table de points
Cette méthode est basée la notion de
Différences divisées, où le polynôme d'interpolation
est obtenu comme ça :
P~(x) =
f [x0] + (x - x0)f [x0, x1] + (x - x0)(x -
x1)f[x0, x1, x2] +
erreur .
Les différences divisées sont obtenues comme
suit :
f[x0] = f(x0)
f[x0] - f[x1]
f [x0,x1] =
~ + (x - x0)(x - x1) +
·
·
·
+ (x - x~_1)f[x0, x1, , x] +
x0 - x1
XO6Xk
Et dans le cas général f[x0,
x1, ... ,xk] = f[xo,xl,...,xk_1]--f[xk_1,xk]
Le programme
//TP4 Calcul Numérique //Polynome De Newton
#include<iostream.h> #include<alloc.h>
#include<conio.h> #include<math.h>
double **NewMat(int n)//Allocation dynamique du matrice
{double * *tmp=(double* *)malloc(2*sizeof(double*));
for(int i=0;i<2;i++)
tmp[i]=(double*)malloc(n*sizeof(double));
return tmp;
}
//Remplissage du tableau des points(Xi,Yi)
void remplir(double **tab,int n)
{ cout<<"Xi:\n";
cout<<"Yi:";
gotoxy(4,wherey()- 1); int x=wherex();
for(int j=0;j<n;j++)
{for(int i=0;i<2;i++)
{cin>>tab[i] [j];
gotoxy(x,wherey());
}
gotoxy(x+=4,wherey()-2);
}
cout<<"\n\n";
}
//Obtention des coefficients du polynome double *coef(double
**tab,int n)
{double *tmp=(double*)malloc(n*sizeof(double)); tmp[0]=tab[ 1]
[0] ;//premier coef
//deuxime coef
tmp[1]=double(tab[1][0]-tab[1][1])/double(tab[0][0]-tab[0][1]);
for(int i=2;i<n;i++)
{double
t=double(tab[1][i-1]-tab[1][i])/double(tab[0][i-1]-tab[0][i]);
tmp[i]=double(tmp[i- 1] -t)/double(tab[0] [0] -tab[0] [i]);
}
return tmp;
}
//Affichage du polynome
void EcrirePoly(double **tab,double *coef,int n) {
cout<<coef[0]<<" ";
for(int i=1 ;i<n;i++)
{if(coef[i] !=0)
{
(coef[i]<0)?cout<<"+("<<coef[i]<<")":cout<<"+"<<coef[i]<<"
"; for(int k=0;k<i;k++) cout<<"(X-"<<tab[0] [k]<<")";
}
}
}
void main()
{clrscr();
cout<<"\t\t\t\tPolynome De Newton\n"; cout<<"Entrer
le nombre des points n ="; int n;
cin>>n;//nombre des points
double **Tabl=NewMat(n);//table des points cout<<"Enter la
table des points\n"; remplir(Tabl,n);
double *cf=coef(Tabl,n);//table des coef cout<<"Le Polynome
de Newton\nP(x)="; EcrirePoly(Tabl,cf,n);
free(Tabl);free(cf);//liberation memoire getch();
}
II.2 - Méthode De Lagrange
Obiectif : Obtention d'un polynôme à
partir d'une taffe de points
_ _
Le polynôme d'interpolation de Lagrange est obtenu par
cette
combinaison :
n
pn (x) = 1 Lk (x)f (xk )
k=0
Où
n
Lk(x) = e x -- xi
xk--xi
i=0,i*k
Le programme
//TP5 Calcul Numerique //Polynome De Lagrange
#include<iostream.h> #include<conio.h>
#include<alloc.h>
double **DefPt(int n)//allocation dynamique de la table des
points
{double **tmp=(double**)malloc(2*sizeof(double*));
for(int i=0;i<2;i++)
tmp[i]=(double*)malloc(n*sizeof(double));
return tmp;
}
//remplissage de la table des points
void remplir(double **tab,int n)
{cout<<"Xi: \n";
cout<<"Yi:";
gotoxy(4,wherey()-1); int x=wherex();
for(int j=0;j<n;j++) {for(int i=0;i<2;i++)
{cin>>tab[i][j];
gotoxy(x,wherey());
} gotoxy(x+=4,wherey()-2);
}
cout<<"\n\n";
}
//calcul et affichage des Lk(x)
int Lk(double **point,int i,int n)
{cout<<"L"<<i<<"(x)=";
double tmp=1;
for(int j=0;j<n;j++)
if(j !=i)
{ cout<<"(x";
if(point[0] [j]<0)
cout<<"+"<<(-1 )*point[0] [j];
if(point[0] [j]>0) cout<<"-"<<point[0] [j];
cout<<")";
tmp*=point[0] [i] -point[0] [j];
}
if(tmp>0)
cout<<"/"<<tmp<<endl;
else
if(tmp<0)
cout<<"/("<<tmp<<")"<<endl;
else
{ cout<<"/"<<tmp<<endl;;
cout<<"\ndivision par zero puisque la condition de points
distincts n'est pas
verifiee"<<endl;
return 0;
}
return 1;
}
//affichage du poly d'interpolation de lagrange
void PolyLagrange(double **point,int n)
{for(int i=0;i<n;i++)
{int ok=Lk(point,i,n);
if(!ok) return;
}
cout<<"Pn(x)=";
for(i=0;i<n;i++)
if(point[ 1] [i])
if(i!=n-1)
cout<<point[ 1] [i]<<"L"<<i<<"(x)+";
else
cout<<point[ 1] [i]<<"L"<<i<<"(x)";
}
void main()
{clrscr();
cout<<"\t\t\t\tPolynome De Lagrange\n";
cout<<"Entrer le nombre des points n ="; int n;
cin>>n;//nombre des points
double **Point=DefPt(n);//table des points cout<<"Enter la
table des points\n";
remplir(Point,n);
cout<<"Le polynome de Lagrange \n";
PolyLagrange(Point,n);
getch();
}
III - Résolution D'une
Equation Non Linéaire De Type f(x) = 0 Dans la pratique, la
plupart des problèmes se ramènent à la
résolution
d'une équation de la forme : f(x) = 0 ,
donc le problème consiste à trouver
la racine dont on sait l'existence et dont, parfois, on
connaît une valeur approchée.
Les méthodes de résolution sont toujours des
méthodes itératives.
III.1 - Méthode De Newton - Raphson Objectif :
résolution de f(x) = 0
Cette méthode s'applique à des équations
de type f(x) = 0, pour
lesquelles on peut calculer la dérivée de f :
f'(x).
Soit x0 une valeur approchée de la racine X inconnue
(x0 doit vérifier la
condition de Fourier f(x0)fgg(x0)
h0).
On a x.1 = x - f(x')
f'(x') ; n: 0, Kmax
Si | x.1 - x | m n alors xn solution
approchée.
Le programme
//TP 6 Calcul Numerique
//Methode de Newton - Raphson //fonction choisie polynome de
degre 3
#include<iostream.h> #include<math.h>
#include<conio.h>
//tableau contient les coefficients du polynome de degre 3
typedef double fct[4];
//tableau contient les coefficients de la 1er derivee
typedef double derive[3];
//tableau contient les coefficients de la 2eme derivee
typedef double derive2 [2];
//condition de Fourier permettant le bon choix de X0 :valeur
initiale
void fourier(fct f,derive2 d2,double &X0)
{
do
cin>>X0;
while ((poly(X0,3 ,f)*poly(X0, 1 ,d2))<0);
}
//methode Newton Raphson Xn+1 =Xn-(F(Xn)/F'(Xn))
//dans ce programme F(x)=fct f et F'(x) = derive d
int N_R(fct f,derive d,double &X0,int Kmax=9,double eps=0.3)
{
double Xnplus1 ,Xn=X0;
for(int i=0;i<Kmax;i++)
{
if(poly(Xn,2,d) != 0)
Xnplus 1=Xn-(poly(Xn,3,f)/poly(Xn,2,d));
else return 0;
if(fabs(Xnplus 1 -Xn)<=eps) {X0=Xnplus 1 ;return 1; }
else Xn=Xnplus1;
}
return 0; }
void main() {clrscr();
cout<<"\t\t\tMethode de Newton Raphson\n";
fct f;
derive d;
derive2 d2;
cout<<"\nSaisir les coefficients a partir du plus haut
degre\n\n"; cout<<"Saisir les coefficients du polynome de degre 3\n";
for(int i=3 ;i>=0;i--) {cout<<"A"<<i<<": ";cin>>f[i]
; } cout<<"\nSaisir les coefficient de la 1 ere derivee\n";
for(i=2;i>=0;i--) {cout<<"A"<<i<<": ";cin>>d[i];}
cout<<"\nSaisir les coefficient de la 2eme derivee\n";
for(i=1 ;i>=0;i--) {cout<<"A"<<i<<":
";cin>>d2 [i];} cout<<"Entrer X0:";
double X0; fourier(f,d2,X0);
int ok;
ok=N_R(f,d,X0);
if(ok) cout<<"Solution ="<<X0;
else cout<<"pas de solution";
getch();
}
III.2 - Accélérateur D 'Aïtken
Objectif : Accélérer la convergence vers
la solution dans le cas de Newton - Rap hson
Dans cette méthode on va essayer
d'accélérer la convergence vers la solution.
Soit x0 donnée, x1=f(x0), x2=f(x1) et par la
suite
xn+1 = xn-2
(xn-1-- n-2)2 X
xn-2 + xn -- 2xn-1
Donc on va construit une suite et à chaque fois on
calcule l'erreur jusqu'on arrive à une erreur < î, mais cette
fois d'une manière plus rapide.
Le programme
//TP7 Calcul Numerique
//Accelerateur D'Aitken
//La fonction choisie exponentielle
#include<iostream.h> #include<alloc.h>
#include<conio.h> #include<math.h> #include<iomanip.h>
typedef double point[3];
//lecture d'un point qui satisfait la condiotion de fourier
void fourier(double &X)
{do
cin>>X;
while((exp(X)*exp(X))<0);
}
void init(point p)//calculer les trois premiers points
{cout<<"Entrer X0 :";
fourier(p[0]);
p[1]=exp(p[0]); p[2]=exp(p[1]); }
double suc(point p)//calculer le prochain point
{double x;
//le calcul avec le coef d'aitken
x=p [0] -(((p[ 1]-p [0])*(p [ 1]-p[0]))/(p [0]+p [2] -(2*p [
1])));
//decalage
p[0]=p[1];
p[1]=p[2];
p[2]=x;
return fabs(p[2]-p[1]);//retourner l'erreur
}
//recherche de la solution
int ait(point p,double eps=0. 1 ,int Kmax=6)
{int k=1;
cout<<"Iteartion Xn-2 Xn-1 Xn L'erreur\n";
cout<<"Iteration 0"<<":";
cout<<setw(6)<<setprecision(3)<<setfill('
')<<p[0]<<" "<<p[1 ]<<" "<<p[2]<<" ";
cout<<fabs(p[2] -p[1 ])<<endl;
if(fabs(p[2]-p[1])<= eps) return 1;
while(k<=Kmax)
{ cout<<"Iteration "<<k<<":";
double tmp=suc(p);
cout<<setw(6)<<setprecision(3)<<setfill('
')<<p[0]<<" "<<p[1 ]<<" "<<p[2]<<" ";
cout<<tmp<<endl;
if(tmp<= eps) return 1;
k++;
}
return 0;
}
void main()
{clrscr();
cout<<"\t\t Acceleration de la convergence"<<endl;
cout<<"\t\t\tMethode d'Aitken"<<endl;
point p;
init(p);//initialisation des 3 points
int ok=ait(p);//s'il ya solution la fonction "ait" retourne 1
if(ok) cout<<"La Solution = "<<p[2];
else cout<<"Pas de solution";
getch(); }
Mais pratiquement on a :
IV - Analyse De Données
Expérimentales
- Méthode Linéaire Des Moindres
Carrés
Objectif : Analyse d'un échantillon de
données relevé expérimentalement et obtenir
l'espérance mathématique de ces données .
Considérons un phénomène physique qui
lie n grandeurs observées xi à y selon une loi linéaire
:
Y = aixi + a2x2 +
·
·
· + anxn
On veut déterminer les coefficients ai par N mesures.
On a théoriquement:
Yi = aixii + a2x2i +
·
·
· + anxni
{
Y2 = aixi2 + a2x22 +
·
·
· + anxn2
·Ym=aixim+a2x2m+
·
·
·
+anxn
m
·
Yi + ri = aixii + a2x2i +
·
·
· + anxni
{
Ym+rm = aixim + a2x2m +
·
·
·+
anxnm
Y2 + r2 = aixi2 + a2x22 +
·
·
· + anxn2
Où les ri sont des erreurs de mesure (résidu) ;
On veut déterminer les ai de sorte que |r| le vecteur de résidus,
soit minimum.
Où 'r' = ri2 +
r22 +
·
·
· +
rm2
Le programme
//TP8 Calcul Numérique //Méthode des Moindres
Carres
#include<iostream.h> #include<iomanip.h>
#include<conio.h> #include<math.h>
#include<alloc.h>
struct PosMax //enregistrement pour le maximum
{int IL;int JC;}; //enregistrement pour la solution
struct solution {double x;int indice;};
////////////////////////////////////////////////////////////////////////////
double **def(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;
}
//calcul du transpose d'une matrice
double **Trans(double **M,int n,int m)
{double * *Mt=def(m,n); for(int i=0;i<m;i++)
for(int j=0;j<n;j++)
Mt[i] [j]=M[j] [i];
return Mt;
}
//produit de 2 matrices
double **Produit(double **M,double **N,int n,int m,int p)
{double * *Q=def(n,p); for(int i=0;i<n;i++)
for(int j=0;j<p;j++)
{double som=0;
for(int h=0;h<m;h++)
som+=M[i] [h] *N[h] [j];
Q [i] [j]=som; }
return Q;
}
//remplissage des echantillons
void remplir(double **M,int n,int m)
{int x=0,y=wherey(); for(int i=0;i<n;i++)
{for(int j=0;j<m;j++)
{gotoxy(x+=8,y);
cin>>M[i] [j]; }
y++;
x=0;
}
}
void print(double **mat,int n,int m)//afficher une matrice
{for(int i=0;i<n;i++) {for(int j=0;j<m;j++)
cout<<setfill('
')<<setw(4)<<setprecision(2)<<mat[i] [j]<<"\t";
//printf("\n");
cout<<endl;
}
}
//////////////////////////Methode de
Gauss////////////////////////////////
double **NewMat(int n,int m) //creation dynamique d'une
matrice
{double * *mat=(double* *)malloc(n*sizeof(double*));
for(int i=0;i<n;i++)
mat[i]=(double*)malloc(m*sizeof(double));
return mat;
}
void charger(double **mat,int n,int m) //initialisation du
matrice
{ cout<<"Entrer la matrice augmentée\n";
for(int i=0;i<n;i++)
{int x=wherex();
for(int j=0;j<m;j++)
{cin>>mat[i] [j] ;gotoxy(x+=6,wherey()- 1); }
cout<<endl;
x=8;
}
}
//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; }
}
}
////////////////////////////////FIN////////////////////////////////////////
//calcul de l'erreur quadratiaue
void residu(double **mat,int n,int m,solution *s,double *r)
{for(int i=0;i<n;i++)
{double som=0;
for(int j=0;j<n;j++) som+=mat[i] [j] *s[j] .x;
r[i]=fabs(som-mat[i] [m-1]);
}
}
void main()
{clrscr();
cout<<"\t\tMethode Des Moindres Carres"<<endl;
cout<<"Entrer le nombre de variables reliee avec
Y"<<endl; int n;cin>>n;
cout<<"Votre Systeme :\nY=C0 X0";
for(int i=1 ;i<n;i++)
cout<<"+C"<<i<<" X"<<i;
cout<<"\nEntrer le nombre d'echantillon :";
int m ;cin>>m;
double **echant=def(m,n+1);
cout<<" Y";
for(i=0;i<n;i++) cout<<" X"<<i;
cout<<endl;
remplir(echant,m,n+1 );//chargement des echatillons
double **M=def(m,n);//definition de M qui contient les valeurs de
X (matrice non carree) for(i=0;i<m;i++) for(int j=0;j<n;j++) M[i]
[j]=echant[i] [j+1] ;//chargement de M
double **Y=def(m,1);//definition du vecteur Y
for(i=0;i<m;i++) Y[i] [0]=echant[i] [0] ;//chargement de Y
double * *Q=Produit(Trans(M,m,n),M,n,m,n);//Q=Mt*M double * *
S=Produit(Trans(M,m,n),Y,n,m, 1 );//S=Yt*Y
//////////////////////Methode de
Gauss////////////////////////////////////
cout<<"\t\t Resolution Du Systeme Linéaire
AX=B\n";
int dim=n;
solution *X=(solution*)malloc(dim*sizeof(solution));
for(i=0;i<dim;i++) X[i] .indice=i;//initialisation des
indices
cout<<endl;
//Declaration de la matrice augmentée double
**A=NewMat(dim,dim+1); //initialisation de A
for(i=0;i<n;i++)
for(j=0;j<n;j++)
A[i] [j]=Q[i] [j];
for(i=0;i<n;i++) A[i] [n]=S [i] [0]; //affichage de A
print(A,n,n+1);
int ok = triangul(A,dim,dim+1 ,X);//triangularisation
cout<<endl;
if(ok)
{ cout<<"Matrice Triangulée\n"; print(A,dim,dim+1
);//A triangulée cout<<endl;
solve(A,dim,dim+1 ,X);//resolution ordre(X,dim);//ordonancement
des solutions
cout<<"Les Solutions Sont :\n"; cout<<endl;
for(i=0;i<dim;i++) cout<<"C"<<X[i]
.indice<<"="<<X[i] .x<<endl;
////////////////////////////////FIN////////////////////////////////////////
double *r=(double*)malloc(dim*sizeof(double));
residu(A,dim,dim+1 ,X,r);//calcul de l'erreur quadratiaue
double R=0;
for(i=0;i<dim;i++) R+=r[i] *r[i];
cout<<"L'erreur Quadratique"<<endl<<"r*r
="<<R;
free(r);
}
free(A);free(X);//liberation memoire
free(Y); free(Q); free(M); free(echant);//liberation memoire
getch();
}
V - Calcul Approché Des
Eléments Propres (valeurs et vecteurs
propres)
- Méthode De La Puissance
Itérée
Objectif : Calcul des valeurs propres sans utiliser le
polynôme caractéristique
La méthode de la puissance itérée est
une méthode itérative de calcul de la valeur propre dominante
d'une matrice (i.e. de plus grand module) et du vecteur propre correspondant.
Soit A une matrice réelle n×n. On supposera que A est
diagonalisable, on notera ëi, (i = 1, n) ses valeurs propres et
Xi, (i = 1, n) ses vecteurs propres correspondant. Les
valeurs propres seront supposées ordonnées:
|ë1| > |ë2| > > |ën|
On a
p = lim
ktu v |lAki|
||~k6~||x
y = lim
k--oo v k
||Ak||)
Ensuite en remplace A par la matrice semblable
B=A-ëXXt, et on reprend le même processus n fois pour
calculer les n éléments propres.
Remarque : on prend k = n+1.
Le programme
//TP9 Calcul Numérique
//Calcul des elements propres //Méthode de La Puissance
Itérée
#include<iostream.h> #include<iomanip.h>
#include<conio.h> #include<math.h> #include<alloc.h>
//structure d'un element propre
struct EP{double val_p;
double *vec_p;};
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(' ')<<A[i] [j];
}
x=0;
y++;
}
}
//calcul de la norme d'un vecteur
double norme(double *v,int n)
{double som=0;
for(int i=0;i<n;i++)
som+=(v[i] *v[i]); return sqrt(som);
}
//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;
}
//calcul de lambda*X*Xt
double **LambdaXXt(double *X,double L,int n)
{double * *tmp=DefMat(n,n);
for(int i=0;i<n;i++) for(int j=0;j<n;j++)
tmp[i] [j]=X[i] *X[j] *L;
return tmp;
}
//calcul de la mtrice semblable B
double **sub(double **A,double **LXXt,int n)
{double * *B=DefMat(n,n);
for(int i=0;i<n;i++)
for(int j=0;j<n;j++)
B[i] [j]=A[i] [j]-LXXt[i] [j];
return B;
}
//calcul des elements propres
EP *ep(double **A,double *V,int n)
{EP *tmp=(EP*)malloc(n*sizeof(EP));
for(int i=0;i<n;i++)
{double *Apmoins 1 =puissance(A,V,n,n),//A puissance P-1
*Ap=puissance(A,Apmoins 1 ,n, 1 ),//A puissance P
NAp=norme(Ap,n),//norme de A puissance P NApmoins 1 =norme(Apmoins 1
,n);//norme de A puissance P-1
//pour eviter la division par zero
if((NAp==0)| | (NApmoins 1==0))
{free(tmp);
return NULL;
}
tmp[i] .val_p= NAp/NApmoins 1 ;//calcul de la valeur propre
//test de distinction des valeurs propres
//ValProp i doit < ValProp i-1
//calcul de vecteur propre
if((i>0)&&(tmp[i] .val_p>=tmp[i- 1] .val_p)) return
NULL; tmp[i] .vec_p=(double*)malloc(n*sizeof(double));
for(int t=0;t<n;t++)
tmp[i] .vec_p[t]=(1 .0/NAp)*Ap[t];
//calcul de la matrice semblable
double **B=A;
A=sub(A,LambdaXXt(tmp[i] .vec_p,tmp[i] .val_p, n),n);
free(B);free(Ap);free(Apmoins 1);
}
return tmp;
}
void main()
{clrscr();
cout<<"\t\tMethode De La Puissance Iteree\n";
cout<<"Donner Le Dimension De La Matrice :";
int n;cin>>n;//dimension de la matrice
double * *A=DefMat(n,n);
cout<<"Charger La Matrice \n";
remplir(A,n,n);
double *V=(double*)malloc(n*sizeof(double));
cout<<"Donner Le Vecteur V \n";
for(int j=0;j<n;j++) cin>>V[j];
EP *ElemProp=ep(A,V,n);
if (ElemProp==NULL) {cout<<"\nImpossible de calculer les
elements propres avec cette methode! ! !\n";
getch();
free(A);free(V);free(ElemProp);
return;}
//affichage des elements propres
for(int i=0;i<n;i++)
{ cout<<"Lambda "<<i<<" =
"<<setprecision(2)<<ElemProp[i] .val_p<<endl; cout<<"Xt
"<<i<<" = "<<"(";
for(j=0;j<n;j++)
cout<<setprecision(2)<<ElemProp[i] .vec_p[j]<<",";
gotoxy(wherex()- 1 ,wherey());
cout<<")"<<endl;
}
free(A);free(V);free(ElemProp); getch();
}
VI - Calcul Approché Des
Coefficients D'un Polynôme de Degré n
Il connu que pour obtenir les valeurs propres d'une matrice
,il faut calculer le déterminant de (A - ë I) ,puis résoudre
le polynôme caractéristique obtenu ;mais si le dimension de la
matrice est grand ,le calcul du déterminant est très fastidieux
.
Pour éviter le problème de calcul du
déterminant on a deux méthodes qui calculent directement les
coefficients du polynôme caractéristique sans passer par le calcul
du déterminant :
- Méthode De Krylov
-Méthode De Sauriau - Faddev
Remarque : dans les deux méthodes
Le premier coefficient (associé au plus haut
degré) est égale à (+1) n Avec la supposition de la
distinction des valeurs propres
|ë1| > |ë2| > > |ën|
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();
}
VI.2 - Méthode De Sauriau - Faddev
Objectif : calcul des coefficients du polynôme
caractéristique
Soit A une matrice carrée; dans cette méthode
les coefficients ai sont obtenus comme suit :
a0= (+1) n
A1=A a1= (+1) n x trace (A1)
B1=A1+a1I
A2=B1A a2= ((+1) n /2) x trace (A2)
B2=A2+a2I
. . .
. . .
. . .
An=Bn+1A an= ((+1) n
/n) x trace (An)
Le programme
//TP1 1 Calcul Numérique
//Calcul des coefficients du polynome caractéristique
//Méthode de Sauriau Faadev
#include<iostream.h> #include<iomanip.h>
#include<conio.h> #include<math.h>
#include<alloc.h>
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++;
}
}
double Trace(double **A,int n)//Calcul de la trace d'une
matrice
{double trace=0;
for(int i=0;i<n;i++) trace+=A[i] [i];
return trace;
}
double **ProduitMat(double **B,double **A,int n,int p,int m)
//Produit matriciel de 2 matrices {double * *tmp=DefMat(n,m);
for(int i=0;i<n;i++)
for(int j=0;j<m;j++)
{double t=0;
for(int k=0;k<p;k++) t+=B[i][k]*A[k][j];
tmp[i] [j]=t;
}
return tmp;
}
void copy(double **A,double **B,int n,int m) //copier matrice A
dans B
{for(int i=0;i<n;i++)
for(int j=0;j<m;j++)
B[i][j]=A[i][j];
}
double **AplustI(double **A,int n,double t) //calcul de B=A+tI
{double * *B=DefMat(n,n);
copy(A,B,n,n);
for(int i=0;i<n;i++) B[i][i]+=t;
return B;
}
int MoinsUn(int dim)//-1 pour degre impaire 1 pour degre paire
{int tmp=1;
if(dim%2) tmp=-1;
return tmp;
}
//calcul des coefficients par methode de Sauriau
double *Sauriau(double **M,int n)
{double *coef=(double*)malloc((n+1 )*sizeof(double));
double * *A=DefMat(n,n);
copy(M,A,n,n);
coef[0]=1 ;//le 1er coefficient
for(int k=1 ;k<=n;k++)
{double tr=-1 *Trace(A,n);
coef[k]=tr/k;
double * *B=AplustI(A,n,coef[k]);
double * *tmp=ProduitMat(B,M,n,n,n);
copy(tmp,A,n,n);
free(tmp);
free(B);
}
free(A);
//produit des coeficients par -1 si degre impaire
// et par 1 si degre paire
int u=MoinsUn(n); if(u==- 1)
for(int i=0;i<n+1 ;i++) coef[i] *=u;
return coef;
}
void main()
{clrscr();
cout<<"\t\t\tMethode de Sauriau Faadev\n\t pour le calcul
des coefficients du Polynome Caracteristique"<<endl;
cout<<"Entrer La Dimension De La Matrice :";int
dim;cin>>dim;
double * *Mat=DefMat(dim,dim);
remplir(Mat,dim,dim);
double *Coef=Sauriau(Mat,dim);
for(int i=0;i<=dim;i++) cout<<"a"<<i<<" =
"<<Coef[i]<<endl;
free(Mat);free(Coef); getch();
}
VII - Recherche De La Valeur Approchée De La
Plus Grande Racine D'un Polynôme de
Degré n
Parmi les problèmes de l'analyse numérique, le
problème de résolution de polynôme de degré >
2.
Pour ça plusieurs méthodes sont apparues, parmi
ces méthode on trouve :
-Méthode De G ra effe
-Méthode De Bernoulli
VII.1 - Méthode De Graeffe
Objectif : Recherche d'une racine d'un
polynôme
On considère un polynôme Pn(x) de
degré n et possède n racines distinctes.
n
Pn(x) = > a)xn6)
t=o
On a si les racines sont éloignées (i.
e|r1 | z |r2 | z z |r3 |) on
peut les calculer de manière approchée comme
suit : r) { - (@
ai_1 ; j: 1, 2
Le but de cette méthode est de rendre les racines
distinctes non- éloignés en des racines éloignées
par la multiplication de f(x) par f (-x)
On a P(x) = (x - a1)(x - a2) ...
(x - an)
P(--x) = (~1)n(x + a1)(x + a2) ... (x +
an)
Alors
P(x) P(--x) = (_1)n(x2 - a1
2)(x2 - a22) ... (x2 -
an2)
Après K itérations on obtient un nouveau
polynôme
n
Qn(y) = > ~)yn6)
t=o
b@
ses racines sont y) { - ; i: 1, 2
b@c5
alors les racines de P(x) sont : x) { }y)
~b ; i: 1, 2
Le programme
//TP12 Calcul Numerique //Methode de Graeffe
#include<iostream.h>
#include<alloc.h>
#include<math.h>
#include<conio.h>
#define carre(x) (x * x)
//le polynome est represente par un tableau
double *DefPoly(int dgr)//allocation dynamique d'un tableau
{double *tmp=(double*)malloc((dgr+1 )*sizeof(double));
return tmp;
}
void LirePoly(double *p,int dgr)//lecture des coefficients du
polynome {for(int i=0;i<=dgr;i++)
{ cout<<"Entrer Le Coefficient a"<<i<<":";
cin>>p[i];
}
}
void PrintSol(double *s,int dgr)
{for(int i=0;i<dgr;i++)
cout<<"La Solution X"<<i<<":
"<<s[i]<<endl;
}
//-1 a la puissance de k
int MoinUn(int k)
{if(k%2) return -1;
else return 1;
}
//transformation du polynome P(x) vers P(x)*P(-x)
void Trans(double *poly,int dgr)
{double *b=DefPoly(dgr);
b[0]=carre(poly[0]);
for(int i=2;i<=dgr+1 ;i++)
{double som=0;
for(int j=1 ;j<=i-1 ;j++)
{if((i+j)>(dgr+1)) break;
som+=MoinUn(j)*poly[i+j-1] *poly[i-j -1];
}
b[i- 1 ]=MoinUn(i- 1 )*(carre(poly[i-1 ])+2*som);
}
for(i=0;i<=dgr;i++) poly[i]=b[i];
free(b);
}
//methode de Graeffe k fois
void graeffe(double *poly,int dgr,int Kmax)
{for(int k=0;k<Kmax;k++) Trans(poly,dgr); }
/*calcul des racines a partir de relation de newton en cas de
coefficients eloignes telque racine Xi=Ai- 1 /Ai /Ai=coef i=0,n */ double
*RelNewton(double *poly,int dgr)
{double *sol=DefPoly(dgr-1);
for(int i=0;i<dgr;i++)
sol[i]=(- 1 )*(poly[ i+1 ]/poly[i]); return sol;
}
//la racine 2*k des solutions
void racine(double *sol,int dgr,int k) {for(int
i=0;i<dgr;i++)
sol[i]=pow(fabs(sol[i]), 1/(2*k));
}
//test de solution
int test(double *p,int dgr,double *sol,double eps=0.1)
{int ok=0;
for(int i=0;i<dgr;i++)
if(fabs(poly(sol[i] ,dgr,p))<=eps)
{ cout<<"Solution = "<<sol[i]<<endl;
return (ok=1);
}
for(i=0;i<dgr;i++)
if(fabs(-poly(sol[i] ,dgr,p))<=eps)
{ cout<<"Solution = "<<-sol[i]<<endl;
return (ok=1);
}
return ok;
}
void main()
{clrscr();
cout<<"\t\tMethode de Graeffe\n";
cout<<"Entrer Le Degre Du Polynome :";int
dgr;cin>>dgr;
double *poly=DefPoly(dgr);
LirePoly(poly,dgr);//lecture des coefs
double *p=DefPoly(dgr);//affecter poly a p pour l'utiliser dans
le test
for(int i=0;i<=dgr;i++) p[i]=poly[dgr-i];
//tantque Kmax augmente le risque de debordement augmente
int ok=0,Kmax=4;
i=1;
//si les coef sont eloignes
double *sol=RelNewton(poly,dgr); ok=test(p,dgr,sol);
//sinon(coefs ne sont pas eloignes) while((!ok) &&
(i<Kmax))
{graeffe(poly,dgr,i);
sol=RelNewton(poly,dgr);
racine(sol,dgr,i);
ok=test(p,dgr,sol);
free(sol);
i++;
}
free(sol);
if(!ok) cout<<"Aucune Solution Trouvee! ! !";
getch();
}
VII.2 - Méthode De Bernoulli
Objectif : Recherche d'une racine d'un
polynôme
Cette méthode est utilisée pour trouver une
racine d'un polynôme
Pn.(x) = aoxn +
aixn +
·
·
· + an
Qui sera transformé en une équation
récurrente aoy(m + n) + aiy(m + n --
1) +
·
·
· + any(m)
Elle est basée sur le principe suivant : on donne les
valeurs suivantes y(0), y(1),..., y (n-1) (solution particulière) et on
va calculer y(n), y (n+1),...
--1
a (aiy(m + n -- 1) +
·
·
·
+ any(m))
o
aoy(m + n) =
On a lorsque m?8, l'expression y (m+1)/y(m) tends vers la
racine du polynôme
Le programme
//TP13 Calcul Numerique //Methode de Bernoulli
#include<iostream.h> #include<conio.h>
#include<alloc.h> #include<math.h>
int bernoulli(double *pol,double *sol,int dgr,double &r,int
Kmax=10,double eps=0.3) {double tmp=0;
double *p=(double*)malloc((dgr+1)*sizeof(double));
//on sauvegarde pol inverse dans p pour utiliser p dans la fct
poly(math.h) for(int t=0;t<=dgr;t++) p[t]=pol[dgr-t];
//on divise les coef par le 1er coef*-1
//afin de construire la suite
for(t=1;t<=dgr;t++) pol[t]*=-(1.0/pol[0]);
//sol[i+1]/sol[i] tend vers la racine
for(t=0;t<dgr;t++)
{if(!sol[t]) continue;
tmp=sol[t+1]/sol[t];
tmp=poly(tmp,dgr,p);
if(fabs(tmp)<eps) {r=sol[t+1 ]/sol[t] ;free(p);return 1; }
}
//calcul des elements de la suite
for(int k=0;k<Kmax;k++)
{tmp=0;
for(int i=1 ;i<=dgr;i++)
tmp+=pol[i] *sol[i- 1];
sol[dgr]=tmp;
tmp=sol[dgr]/sol[dgr- 1];
tmp=poly(tmp,dgr,p);
if(fabs(tmp)<eps) {r=sol[dgr]/sol[dgr- 1] ;free(p);return 1;
}
else
for(i=0;i<dgr;i++)
sol[i]=sol[i+1 ];
}
cout<<"pas de convergence vers la solution";
free(p);
return 0;
}
void main()
{cout<<"Entrer le degre du polynome n :";
int dgr;cin>>dgr;//degre du poly
cout<<"Entrer les coefficient a partir du plus haut
degre\n";
double *pol=(double*)malloc((dgr+1 )*sizeof(double));
//lecture des coef
for(int i=0;i<=dgr;i++)
{ cout<<"a"<<i<<"=";
cin>>pol[i];
}
double * SolPart=(double*)malloc((dgr+1 )*sizeof(double));
cout<<"Entrer la solution particuliere\n";
//lecture de la solution particuliere
for(i=0;i<dgr;i++)
{ cout<<"s"<<i<<"=";
cin>>SolPart[i];
}
double racine;
if(bernoulli(pol,SolPart,dgr,racine)) cout<<"racine
approchee ="<<racine; free(pol);free(SolPart);//liberation memoire
getch();
}
| 
