Projet de 2 semaines du 27 janvier au 4 fevrier

I. Introduction *

II. Algorithme de PA=LU *

III. Cas des matrices creuses *

IV. Cas des matrices tridiagonales *

V. Le problème des approximations *

VI. Que fait notre programme ? *

VII. Script de notre algorithme ( au format MatLab ) *

1. Introduction
1. Notre choix

Après réflexion, nous avons choisi le sujet numéro 1 qui traite le problème de la factorisation d’une matrice A par le calcul de PA=LU .

Nous utiliserons comme logiciel de travail MatLab© .Tous nos algorithmes qui vont être décrit par la suite pourront être visionnés sur ce logiciel 

2. Le SUJET 1

La méthode d'inversion de systèmes et de matrices dite méthode d'élimination de Gauss peut s'écrire pour des problèmes carrés comme une factorisation de la matrice A du problème sous la forme P A = L U où L et U sont des matrices triangulaires inférieure et supérieure ( L ayant une diagonale formée d'éléments unité) et P est une matrice de permutation.

 

Vous écrirez un algorithme qui permet de réaliser cette factorisation, et de l'appliquer à la résolution de systèmes linéaires et l'inversion de la matrice A. Cet algorithme sera porté sur le logiciel MATLAB.

Vous compterez (analytiquement et pratiquement) le nombre d'opérations - additions et multiplications - nécessaire à la factorisation et à l'inversion, en essayant de le minimiser, et observerez l'effet des erreurs d'arrondi sur la précision des résultats.

Enfin, vous considérerez le cas des matrices creuses ( ayant un grand nombre de 0 ) et en particulier celui des matrices tridiagonales. Ici encore, vous essaierez de tenir compte de la structure creuse des matrices afin de réduire le nombre d'opérations ainsi que la place mémoire nécessaire.



2. Algorithme de PA=LU
1. Algorithme utilisé

Déclarations

Soient

A une matrice donnée d’ordre N (créée de manière aléatoire)

P matrice des permutations (Ps (ei)=es (i)) ( inconnue )

L matrice triangulaire inférieure ayant des 1 sur la diagonale ( inconnue )

U matrice triangulaire supérieure ( inconnue )

C ,E sont des matrices temporaires

P,L,A sont initialisés à la matrice identité

P, L et U sont de la même dimension que A.

On commence l’algorithme à partir de la première ligne ( => ligne courant i est égale à 1 )

Soit A’ la sous matrice de A définie comme ci-après et qui dépend de la ligne en cours

A =

a b c d

e f g h

i j k l

m n o p

A’ = (si deuxième ligne en cours)

f g h

j k l

n o p

 

Tant que la sous matrice A’ n’est pas triangulaire inférieure

Première étape

On initialise L et C (matrice temporaire) à l’identité

On cherche d’abord un pivot qui soit diffèrent de zéro

Pour un choix convenable de ce pivot , ce référer au chapitre `Le problème de approximations`

Si le pivot se trouve sur une autre ligne p ( i<>p)

Dans L on fait : Ligne(i) <-> Ligne(p)

Dans A on fait : Ligne(i) <-> Ligne(p)

Deuxième étape

On calcule P = P * L ( ce qui permet de mettre à jour la matrice de permutation des lignes).

Troisième étape

On applique la méthode du pivot de gauss pour faire apparaître des zéros sur la colonne courante i ( sauf le premier élément). Pour cela on utilisera une boucle

Soit j ( i+1 < j < N) le numéro de la ligne à modifier

Si on arrive au cas où L(i)<->L(i)+aL(j) alors

- On calcul le coefficient ( -A(j,i)/A(i,i) ) = a

- On met ce coefficient dans C(j,i)

- On fait la combinaison linéaire :

L(i)<->L(i)+aL(j)

On calcul maintenant E : E = C * L * E

Quatrième étape

On recommence l’étape 1 tant que A n’est pas triangulaire supérieure

On calcul L qui est égal à P*inv(E)

U est égal à A ( car A a été modifier à chaques étapes de l’algorithme et est maintenant triangulaire supérieure )

Et P a été calculé au cours des étapes précédentes

2. But de cet algorithme

Ce que nous voulons faire, c’est mettre une matrice A quelconque ( avec det(A)<>0 entre autre ) sous la forme de produit de matrices triangulaires, ce qui simplifierait les calculs de type résolution d’un système Ax=b et inversion d’une matrice.

Pour cela on va déterminer trois matrices

P matrice de permutation , qui permet de savoir quelles permutations de lignes ont eu lieu.

U matrice triangulaire supérieure

L matrice triangulaire inférieure contenant que des 1 sur la diagonale

 

telles que P * A = L * U

3. Résolution Ax=b

A matrice carrée telle que P A = L U

b matrice colonne N*1

Grâce à la formule PA=LU, on peut écrire

Soit y une matrice colonne

y = L^-1 * P * b

x = U^-1 * y

Remarque : pour l’inversion d’une matrice nous n’utilisons pas la fonction prédéfinie de MatLab . Nous utiliserons une fonction d’inversion qui sera spécifique à la matrice à inverser ( triangulaire on diagonale .Cette fonction sera utilise indirectement pour la résolution de Ax=b et pour l’inversion de A, qui est le point suivant)

4. Inversion de A

A matrice carrée telle que P A = L U ( donc A est inversible)

On a P A = L U

=> P = L U A^-1

=> L^-1 P = U A^-1

=> A^-1 = U^-1 L^-1 P

5. Remarques

- L’inversion de matrices triangulaires étant plus rapide cette méthode de factorisation a un gros avantage

- En effet , pour inverser une matrice quelconque , on peut utiliser la méthode de Gauss-Jordan , qui fait en gros deux pivots de gauss ( un en descendant et un en montant )

- Comme U et L sont triangulaires , il n’y aura ici qu’un pivot de Gauss à faire ( montant pour U et descendant pour L)

- La multiplication de matrices triangulaires est aussi facilité car il y a beaucoup de zéros

 

6. Complexité

Lorsque des opérations sont effectués dans une boucle la complexité est multiplié par n et par n² si elle se situe deux boucles imbriquées

Le compteur interne compte toutes les opérations test divers et calcul sur matrice

Un compteur subsidiaire existe il correspond au nombre d’affectation ( " place " mémoire ) nous avons essayé de l’utiliser tout au long de nos essais

 

Pour le calcul de PA=LU

Recherche du pivot max

il y a un parcourt d’une colonne et comparaison de valeurs

si on considère que la ligne courant est i , il y a alors (n-i+1) éléments parcourus et (n-i+1) comparaison .

il y a 3 opérations par éléments à permuter , et ceci pour (n-i) éléments

ce qui fait 2*3(n-i) affectations *2 pour la matrice de P.

Produit de deux matrices carrées

la complexité est de n²(2n-1) ~ 2n³

Pour une ligne de A il y a

(n) multiplications

(n-1) additions.

On a donc (2n-1) opérations par colonnes de B.

Or il y a n colonnes et n lignes ce qui fait nxnx(2n-1) opérations.

Produit d’ une matrice carrée et un vecteur colonne

nx(2n-1) ~ 2n²

Opération du pivot : L(i)<->L(i)+aL(j)

il y a (n-i) additions et (n-i) multiplications

on a aussi (n-i) affectations

on a donc ici 2(n-i) ( additions + multiplications )

Inversion d’une matrice quelconque

pour la transformation en triangulaire supérieure 5 n³

pour la diagonale unitaire 2n

pour la transformation en identité 5 n³

Au total 5n²(2n+1)

Coût total de la recherche de P , L et U

6n⁴+14 n³+11 n²

Pour la résolution Ax=b

Inversion d’une matrice triangulaire

pour la diagonale unitaire 2n

pour la transformation en identité 5 n³

Au total 5n³+2n

Coût total

12 n³ +4 n²+4n

Pour inv(A)

Coût total

invU*invU*P

2*(5n³+2n)+2*n³

Au final

calcul de LU + invA

6n⁴+24n³+o(n²)



3. Cas des matrices creuses
1. Méthode d’optimisation

Notre méthode consiste ici à éliminer tous les calculs inutiles (i.e. quand nous multiplions pas zéro ou on soustrait 0 )

En fait on reprend l’algorithme précèdent avec quelques légères modifications . Dans la boucle de la troisième étape , si le coefficient A(j,i) est égal à zéro ,on ne fait aucun calcul (i.e. L(i) <-> L(i)+aL(j) ).

2. Complexité

Elle dépend du nombre de zéro qu’il y a dans la matrice initiale



4. Cas des matrices tridiagonales
1. Méthode d’optimisation

Notre méthode consiste à utiliser l’algorithme de Hessenberg qui a l’avantage d’être spécialement conçu pour ce type de matrices

Que nous dit cet algorithme ?

Soit A =

b ₁ c ₁

a ₂ b₂ c₂

.............

................

a_n-1 b _n-1 c _n-1

a _n b _n

 

On cherche L et U définit pas

Soit L =

1

l ₂ 1

..........

...........

l_n-1 1

l_n 1

 

Soit U =

d ₁ u ₁

d₂ u₂

........

..........

d _n-1 u _n-1

d _n

 

Analytiquement nous avons les propriétés suivantes:

b₁ = d₁

c_i = u_i (2 <= i <= N)

l_i = a_i / d_i-1 (-1-)

d_i = b_i - l_i * c_i-1

P est ici la matrice identité

Les résultats que nous avons observé

Nous pouvons voir que le temps de calcule de L et U est court En effet nous ne traitons que les coefficients nécessaires et nous ne faisons aucun calculs inutiles

Pour la résolution A x = f

En conséquence de la méthode de Hessenberg et en suivant ces deux étapes :

- L y = f

- U x = y

Nous arrivons aux résultats suivants :

y₁ = f₁ (-2-)

y_i = f_i - l_i * y_i-1 ( 2 <= i <= N )

x_N = y_N / d_N (-3-)

x_i = ( y_i - c_i * x_i+1 ) / d_i ( (N-1) >= i >= 1 )

2. Complexité

Pour le calcul de L et U est de ( N-1)(addition + multiplication + division) . ( trivial : cf. formules (-1-) )

Pour le calcul des y_i ( cf. (-2-) ) nous avons ( N-1)(addition + multiplication)

Pour le calcul des x_i nous avons ( N-1)(addition + multiplication) et N(division) ( cf. (-3-) )

Donc pour la résolution Ax=f il faut

3( N-1)(addition + multiplication) + (2N-1)(division)



5. Le problème des approximations
1. Cas général

Le problème vient du choix des pivots lors des procédures de  " Gaussisation " de matrices

En effet , théoriquement , si nous prenons un pivot trop petit en valeur absolue , le résultat fournit par un pivot de Gauss peut être carrément faux , cette erreur est due aux arrondis faite par les calculatrices et les ordinateurs sur les nombres décimaux. C’est pour cela que nous avons choisi de prendre comme pivot le coefficient le plus grand en valeur absolue.

Dans notre programme nous avons calculé l’inverse de A de deux manières

La première avec le choix d’un pivot maximum pour le pivot

La deuxième avec le choix d’un pivot minimum pour le pivot

Nous avons fait apparaître qu’il y avait quand même des différences sur les résultats finaux

Par contre avec tous nos tests , nous n’avons pas trouvé de cas où nous avions des erreurs supérieures à 10^-14

2. Cas des matrices creuses

En utilisant la méthode de Hessenberg ce problème est beaucoup moins important car ici nous n’utilisons pas de pivot . Donc la précision des calculs est bonne

Par contre en utilisant la méthode du pivot de Gauss nous avons une moins bonne précision, car il faut à chaque fois choisir entre deux pivots.

Il s’en suit un gain de temps , de performances d’exécution , de précision au niveau des calculs. En faite : utiliser un algorithme spécifique à une tache précise pour optimiser...



6. Que fait notre programme ?
1. Utilisation

il faut exécuter le fichier appli.M avec MatLab

on demande à l’utilisateur de donner une dimension N, qui servira à créer et initialiser les diverses matrices à l’ordre N

Il s’agit plus particulièrement des matrices tridiagonales

En effet afin de mieux percevoir les effets d’erreurs d’arrondi nous utilisons comme résultats de référence ceux calculés par la méthodes de Henssenberg

Méthode d’Henssenberg

L’ordinateur calcul les matrices L et U ainsi que les solution de AX=b

Comparaison du nombre d’opérations effectuées par le compteur interne de notre algorithme et le compteur de matlab FLOPS

Méthode du pivot Maxi

Calcul des matrices L , U , P ainsi que les solutions de AX=b

Comparaison du nombre d’opérations

Calcul de l’inverse de A

L’inverse est calculée de 3 manières

Par l’inverse de U et L * P

Par la fonction de mathlab

Par la fonction invA " double pivot "

Vérification de PA=LU

Affichage des matrices P*A et L*U et Lh *Uh ( issu de la méthode de Henssenberg )

Affichage des solutions

Affichages des solutions issues des différentes méthodes

Pivot Maximum

Henssenberg

Pivot aléatoire Minimum

Comparaison des solutions

Différences entre les solutions issues de différentes méthodes avec les solutions de la méthode Henssenberg

2. Limites

Utilisation du module inv_trig.m

Spécialement conçu pour inverser des matrices triangulaires

Il optimise le nombre des opérations du fait de la forme particulière des matrices considérées

Utilisation de la transposée pour les matrices triangulaires inférieures

 

3. Résultat à l’écran

 

" appli

Pour une présentation des résultats satisfaisants prendre une matrice <7 Remarque pour les erreurs d arrondi privilégier une matrice plus grande

Entrer la dimension de la matrice :5

L ordinateur va générer une matrice aléatoire de dimension n tridiagonale

Calcul D OPTIMISATION

REMARQUE : il existe une différence notable entre le compteur interne et la fonction FLOPS cela est du essentiellement à la mise en place

du compteur interne en effet la fonction FLOPS compte également

les opérations d incrémentations de ce compteur

Nombre d opérations pour le calcul des solutions par la méthode de Hensenberg avec la fonction FLOPS

flops1 =

224

Nombre d opérations pour le même calcul avec un compteur interne

noph =

33

Nombre d opérations pour le calcul des solutions par la methode du PIVOT

avec la fonction FLOPS

flops2 =

6101

Nombre d opérations pour le même calcul avec un compteur interne

nopp =

4932

Pressez une touche

 

Calcul et vérification de l inverse

INVA_LUP =

-2.8139e+000 1.5787e+000 -5.8853e-001 -2.1879e+000 1.8182e+000

3.3874e+000 -1.3362e-001 4.9813e-002 1.8518e-001 -1.5389e-001

-5.4929e-001 2.1667e-002 4.1222e-001 1.5324e+000 -1.2735e+000

-1.2142e+000 4.7895e-002 9.1121e-001 -2.3138e+000 1.9229e+000

3.3543e+000 -1.3231e-001 -2.5173e+000 6.3919e+000 2.1596e-001

inv_matlab =

-2.8139e+000 1.5787e+000 -5.8853e-001 -2.1879e+000 1.8182e+000

3.3874e+000 -1.3362e-001 4.9813e-002 1.8518e-001 -1.5389e-001

-5.4929e-001 2.1667e-002 4.1222e-001 1.5324e+000 -1.2735e+000

-1.2142e+000 4.7895e-002 9.1121e-001 -2.3138e+000 1.9229e+000

3.3543e+000 -1.3231e-001 -2.5173e+000 6.3919e+000 2.1596e-001

inv_fct =

-2.8139e+000 1.5787e+000 -5.8853e-001 -2.1879e+000 1.8182e+000

3.3874e+000 -1.3362e-001 4.9813e-002 1.8518e-001 -1.5389e-001

-5.4929e-001 2.1667e-002 4.1222e-001 1.5324e+000 -1.2735e+000

-1.2142e+000 4.7895e-002 9.1121e-001 -2.3138e+000 1.9229e+000

3.3543e+000 -1.3231e-001 -2.5173e+000 6.3919e+000 2.1596e-001

Pressez une touche

 

VERIFICATION DE PA=LU

P*A est egale à Lh*Uh calculer par Henssenberg aux permutations de lignes près !! Probablement dues aux permutations lors de la recherche du pivot max

LhUh =

2.6876e-002 3.1754e-001 0 0 0

6.8135e-001 7.0982e-001 8.8699e-001 0 0

0 3.8581e-001 9.3790e-001 6.5206e-001 0

0 0 3.8773e-001 2.3991e-001 1.5034e-001

0 0 0 4.9974e-001 1.8090e-001

 

 

LU est bien egale PA par la décomposition du pivot

PA =

6.8135e-001 7.0982e-001 8.8699e-001 0 0

0 3.8581e-001 9.3790e-001 6.5206e-001 0

2.6876e-002 3.1754e-001 0 0 0

0 0 0 4.9974e-001 1.8090e-001

0 0 3.8773e-001 2.3991e-001 1.5034e-001

LU =

6.8135e-001 7.0982e-001 8.8699e-001 0 0

0 3.8581e-001 9.3790e-001 6.5206e-001 0

2.6876e-002 3.1754e-001 -5.6379e-018 2.1684e-019 0

0 0 0 4.9974e-001 1.8090e-001

0 0 3.8773e-001 2.3991e-001 1.5034e-001

Pressez une touche

Solution de AX=b par le pivot MAX

x1 =

3.617583378963669e-001

3.801244433953384e-001

-4.795001040808731e-001

8.858902018682369e-001

-1.571203076753021e-001

Solution de AX=b par Henssenberg

x =

3.617583378963673e-001

3.801244433953384e-001

-4.795001040808729e-001

8.858902018682366e-001

-1.571203076753021e-001

Solution de AX=b par le pivot MIN

x2 =

3.617583378963458e-001

3.801244433953404e-001

-4.795001040808739e-001

8.858902018682369e-001

-1.571203076753022e-001

 

Calcul des erreurs d arrondi prenons comme base les solutions calculées par Henssenberg différence x-x1 c-a-d Hens-Pivot MAX

ans =

3.885780586188048e-016

0

1.665334536937735e-016

-2.220446049250313e-016

2.775557561562891e-017

différence x-x2 c-a-d Hens-Pivot MIN

ans =

2.142730437526552e-014

-1.998401444325282e-015

9.992007221626408e-016

-2.220446049250313e-016

8.326672684688673e-017

CONCLUSION la plus petite différence est celle entre les valeurs

de référence et les valeurs calculées par le pivot Max en effet lors des calculs on multiplie par un coef=element pivot/élément à supprimer plus l élément pivot est grand et plus le coef multiplicateur est précis

"



7. Script de notre algorithme ( au format MatLab )
1. Script du module principal : appli.m

clear;

clc;

flops(0);

format short e;

nop=0;

naf=0;

[n,nop,naf,A,b]=init(nop,naf);

 

if (det(A)~=0)

[x,Lh,Uh,noph,nafh]=lu_hens(A,b,n,nop,naf);

flops1=flops;

disp(' ')

disp('Calcul D OPTIMISATION')

disp('REMARQUE : il existe une difference notable entre le compteur interne et')

disp('la fonction FLOPS cela est du essentiellement à la mise en place ')

disp('du compteur interne en effet la fonction FLOPS compte également ')

disp('les opérations d incrémentations de ce compteur ')

disp('Nombre d opérations pour le calcul des solutions par la methode de')

disp('Hensenberg avec la fonction FLOPS ')

flops1

disp('Nombre d opérations pour le meme calcul avec un compteur interne ')

noph

[x1,L,U,P,nopp,nafp]=lu_pivom(A,b,n,nop,naf,1);

disp('Nombre d opérations pour le calcul des solutions par la methode du PIVOT')

disp(' avec la fonction FLOPS ')

flops2=flops-flops1

disp('Nombre d opérations pour le meme calcul avec un compteur interne ')

nopp

disp(' Pressez une touche ')

disp(' ')

pause

clc

 

%% INVERSION DE A %%

disp('Calcul et verification de l inverse');

INVA_LUP =inv(U)*inv(L)*P

inv_matlab=inv(A)

[invA,k,l]=inv_A(A,n,0,0);

inv_fct=invA

disp(' Pressez une touche ')

disp(' ')

pause

clc

 

disp('VERIFICATION DE PA=LU')

disp('P*A est egale à Lh*Uh calculer par Henssenberg aux permutations de lignes pres !! ')

disp('Probablement dues aux permutations lors de la recherche du pivot max')

LhUh=Lh*Uh

disp('LU est bien egale PA par la décomposition du pivot')

PA=P*A

LU=L*U

disp(' Pressez une touche ')

disp(' ')

pause

clc

 

format long e;

disp('Solution de AX=b par le pivot MAX')

x1

disp('Solution de AX=b par Henssenberg')

x

[x2,L,U,P,l,k]=lu_pivom(A,b,n,l,k,0);

disp('Solution de AX=b par le pivot MIN')

x2

pause

clc

 

format long e;

disp('Calcul des erreurs d arrondi prenons comme base les solutions calculées par Henssenberg')

disp('difference x-x1 c-a-d Hens-Pivot MAX');

x-x1

disp('difference x-x2 c-a-d Hens-Pivot MIN');

x-x2

disp('CONCLUSION la plus petite difference est celle entre les valeurs de reference')

disp('et les valeurs calculées par le pivot Max en effet lors des calculs on multiplie ')

disp('par un coef=element pivot/element a suprimer plus l element pivot est grand et ')

disp('plus le coef multiplicateur est precis ')

 

 

else

disp('Resolution Impossible')

end

 

2. Script du module d’initialisation : init.m

function [n,nop,naf,A,b]=init(nop,naf)

disp('Pour une presentation des résultats satisfaisants prendre une matrice <7')

disp('Remarque pour les erreurs d arrondi privilégier une matrice plus grande ')

n=input('Entrer la dimension de la matrice :');

disp('L ordinateur va générer une matrice aleatoire de dimension n tridiagonale')

b=rand(n,1);

naf=naf+n;

A=zeros(n);

naf=naf+(n*n);

% INITIALISATINO DE MATRICE A TRIGIAGONALE %

for i=1:n

A(i,i)=rand; % b(i)=A(i,i)

naf=naf+1;

end

for i=1:n-1

A(i,i+1)=rand; % c(i)=A(i,i+1)

naf=naf+1;

end

for i=2:n

A(i,i-1)=rand; % a(i)=A(i,i-1)

naf=naf+1;

end

 

3. Script du module de permutation : per_pivo.m

function [A,P,nop,naf]=per_pivo(A,P,n,i,nop,naf)

 

% i= colonne de recherche

%% RECHERCHE DU PIVOT OPTIMISE DU MAX %%

p=i;

naf=naf+1;

tmp=abs(A(i,i));

naf=naf+1;

nop=nop+2;

for j=i+1 : n

nop=nop+2;

if (abs(A(j,i))>tmp)

p=j;

naf=naf+1;

tmp=abs(A(j,i));

nop=nop+1;

naf=naf+1;

end

end

nop=nop+1;

if (tmp==0)

disp('Resolution impossible')

break;

end

 

%% PERMUTATION ISSUE DU PIVOT MAX %%

Ltmp=A(p,1:n);

naf=naf+n;

A(p,1:n)=A(i,1:n);

naf=naf+n;

A(i,1:n)=Ltmp;

naf=naf+n;

Ltmp=P(p,1:n);

naf=naf+n;

P(p,1:n)=P(i,1:n);

naf=naf+n;

P(i,1:n)=Ltmp;

naf=naf+n;

 

4. Script du module d’inversion :inv_A.m

function [P,nop,naf]=inv_A(A,n,nop,naf)

 

P=eye(n,n);

naf=naf+n;

 

%% pour arriver à A matrice triangulaire supérieure

for i=1 :n-1

[A,P,naf,nop]=per_pivo(A,P,n,i,naf,nop);

nop=nop+1;

if (A(i,i)~=0)

for j=i+1 : n

coef=-A(j,i)/A(i,i);

P(j,1:n)=(coef*P(i,1:n))+P(j,1:n);

A(j,1:n)=(coef*A(i,1:n))+A(j,1:n);

naf=naf+n+n+1;

nop=nop+2+2*n+2*n;

end

end

end

 

%% Matrice A avec diagonale unitaire

for i=1:n

nop=nop+1;

if (A(i,i)~=0)

P(i,1:n)=P(i,1:n)/A(i,i);

nop=nop+n;

naf=naf+n;

A(i,1:n)=A(i,1:n)/A(i,i);

nop=nop+1;

naf=naf+1;

end

end

 

%% pour arriver à A matrice identité

for k=0 :n-2

i=n-k; % i 2 a n

for j=1 : i-1

nop=nop+1;

if (A(i,i)~=0)

coef=-A(j,i)/A(i,i);

naf=naf+1;

nop=nop+1;

P(j,1:n)=(coef*P(i,1:n))+P(j,1:n);

nop=nop+n*2;

naf=naf+n;

A(j,1:n)=(coef*A(i,1:n))+A(j,1:n);

nop=nop+2*n;

naf=naf+n;

end

end

end

 

 

 

5. Script du module LU par pivot : LU_pivom.m

function [x1,L,U,P,nop2,naf2]=lu_pivom(A,b,n,nop,naf,MAX);

% CALCUL DE LA DECOMPOSITION LU=PA AVEC PIVOT MAX

 

% CALCUL DE LA DECOMPOSITION LU=PA AVEC PIVOT MAX

P=eye(n,n);

naf=naf+n;

E=eye(n,n);

naf=naf+n;

X=eye(n,n);

naf=naf+n;

for i=1: n-1

L=X;

naf=naf+(n*n);

C=X;

naf=naf+(n*n);

%% RECHERCHE DU PIVOT OPTIMISE DU MAX %%

p=i;

naf=naf+1;

tmp=abs(A(i,i));

naf=naf+1;

nop=nop+1;

if (MAX==1)

for j=i+1 : n

nop=nop+2;

if (abs(A(j,i))>tmp)

p=j;

naf=naf+1;

tmp=abs(A(j,i));

naf=naf+1;

nop=nop+1;

end

end

else

if (tmp==0)

tmp=2e50 % gestion du premier pivot nul

end

for j=i+1 : n

nop=nop+2;

if (abs(A(j,i))<tmp) & (abs(A(j,i))~=0)

p=j;

naf=naf+1;

tmp=abs(A(j,i));

naf=naf+1;

nop=nop+1;

end

end

end

nop=nop+1;

if (tmp==0) | (tmp==2e50)

disp('Resolution impossible ');

break;

end

%% PERMUTATION ISSUE DU PIVOT MAX %%

Ltmp=A(p,i:n);

naf=naf+(n-i);

A(p,i:n)=A(i,i:n);

naf=naf+(n-i);

A(i,i:n)=Ltmp;

naf=naf+(n-i);

 

Ltmp=L(p,i:n);

naf=naf+(n-i);

L(p,i:n)=L(i,i:n);

naf=naf+(n-i);

L(i,i:n)=Ltmp;

naf=naf+(n-i);

 

P=L*P;

naf=naf+(n*n);

nop=nop+(2*n-1)*n*n;

for j=i+1:n

nop=nop+1;

if (A(j,i)~=0)

C(j,i)=-A(j,i)/A(i,i);

nop=nop+1;

naf=naf+1;

A(j,i:n)=(C(j,i)*A(i,i:n))+A(j,i:n);

naf=naf+(n-i);

nop=nop+2*(n-i);

end

end

E=C*L*E;

naf=naf+(n*n);

nop=nop+(2*n-1)*n*n+(2*n-1)*n*n;

end

 

[invE,nop,naf]=inv_A(E,n,nop,naf);

L=P*invE;

%L=P*inv(E);

naf=naf+(n*n);

nop=nop+(n*n);

U=A;

naf=naf+(n*n);

 

%% RESOLUTION DU SYSTEME %%

[invL,nop,naf] =inv_A(L,n,nop,naf);

Y=invL*P*b;

%Y=inv(L)*P*b;

naf=naf+n;

nop=nop+(2*n-1)*n*n+(2*n-1)*n;

[invU,nop,naf] =inv_A(U,n,nop,naf);

x1=invU*Y;

%x1=inv(U)*Y;

naf=naf+n;

nop=nop+(2*n-1)*n;

nop2=nop;

naf2=naf;

 

6. Script du module factorisation LU par HENSSENBERGl : lu_hens.m

function [x,L,U,nop1,naf1]=LU_HENS(A,b,n,nop,naf)

L=eye(n,n);

naf=naf+n;

U=L;

naf=naf+(n*n);

% CALCUL DES MATRICES L ET U METHODE HESSENBERG %

U(1,1)=A(1,1);

naf=naf+1;

for i=1 : n-1

U(i,i+1)=A(i,i+1); % u(i)=c(i)

naf=naf+1;

end

for i=2:n

L(i,i-1)=A(i,i-1)/U(i-1,i-1);

naf=naf+1;

nop=nop+1;

U(i,i)=A(i,i)-L(i,i-1)*A(i-1,i);

naf=naf+1;

nop=nop+2;

end

 

% RESOLUTION DE AX=b %

% 1er étape Ly=b %

y(1,1)=b(1,1);

nap=nap+1;

for i=2 : n

y(i,1)=b(i,1)-L(i,i-1)*y(i-1,1);

naf=naf+1;

nop=nop+2;

end

 

% 2eme étape Ux=y %

x(n,1)=y(n,1)/U(n,n);

naf=naf+1;

nop=nop+1;

for k=1 : n-1

i=n-k;

x(i,1)=(y(i,1)-A(i,i+1)*x(i+1,1))/U(i,i);

naf=naf+1;

nop=nop+3;

 

end

naf1=naf;

nop1=nop;

7. Script du module d’optimisation : inv_trig.m
1. L’inverse pour les matrices triangulaires ( non utilisé ... )