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|
| 
