TP MAO 5

> restart: with(linalg): with(pldx):

>

Cas où l'algorithme converge

> ma:= randmatrix(4,4);
# les coefficients vont de -99 à +99

> ma := matrix([[-53, 85, 49, 78], [17, 72, -99, -85], [-86, 30, 80, 72], [66, -29, -91, -53]]);

> v0:= evalf(randvector(4));
# on se place sur les réels et pas dans le corps des fractions

> itera:= proc(v1); evalm(ma &* v1): evalm(%/norm(%) ); end;

> (itera@@30)(v0); itera(%); vp1:= norm(evalm(ma &* %));
# la direction ne change plus beaucoup

> Digits:=20: (itera@@190)(v0); itera(%); vp1:= norm(evalm(ma &* %)); Digits:=10:
# on augmente la précision

190 étapes pour gagner 18 décimales, cela fait un taux de convergence égal à :

> (10.^18)^(1/190);

Et cela se vérifie...

> [eigenvals(map(evalf, ma))]; sort(map(abs,%)); %[-1]/%[-2];

Cas où le premier algorithme ne converge pas

> ma:= randmatrix(4,4):

> ma := matrix([[-11, 38, -7, 58], [-94, -68, 14, -35], [-14, -9, -51, -73], [-73, -91, 1, 5]]);

> v0:= evalf(randvector(4));

> itera:= proc(v1) global ma ; evalm(ma &* v1): return evalm(%/norm(%) ); end;

> (itera@@90)(v0); itera(%); vp1:= norm(evalm(ma &* %));
# la direction ne se stabilise pas

Il y a donc plusieurs valeurs propres dominantes.
La matrice étant réelle, il y a des valeurs propres conjuguées.

Diviser par la norme (qui est un réel) ne permet plus d'identifier la direction (colinéarité avec coefficient complexe). On selectionne la composante dominante

> cpodom:= proc(w) local modmx, jx, j, tmp;
modmx:= 0; jx:=1;
for j from 1 to 4 do
tmp:= evalf(abs(w[j])); if tmp > modmx then modmx:=tmp; jx:=j; fi;
od; return w[jx];
end :

On teste

> evalm(randvector(4)+randvector(4)*I); cpodom(%);

Réécriture de l'itérateur

> iterb:= subs(ma=mb, norm=cpodom, eval(itera));

Décalage de la matrice initiale

> tr:=100*I; mb:= evalm(ma+tr);

Convergence : la direction se stabilise

> (iterb@@30)(v0); iterb(%); vp1:= cpodom(evalm(mb &* %));

Augmentation de la précision et du nombre d'étapes

> Digits:=20: vn:= (iterb@@50)(v0); vn1:= iterb(%); vp1:= cpodom(evalm(mb &* %));
evalm(vn1-vn); norm(%,2); Digits:=10:

Vitesse : 11 décimales en 50 étapes

> (10.^11)^(1/50);

Vérification

> [eigenvals(map(evalf, mb))]; sort(map(abs,%)); %[-1]/%[-2];

Les valeurs propres dominantes de ma sont donc

> mu1:= evalf(vp1-tr, 20); mu2:= conjugate(%);

Et les vecteurs propres

> wp1:= pldx[xpurge](vn): Vector[column](wp1);
ma &* wp1 - 'mu1'*wp1; evalf(evalm(%), 20): norm(%,2);
# remarque sur les arrondis de calcul : on perd deux décimales sur le résultat

Calcul de la valeur propre de plus petit module

Réécriture de l'itérateur

> iterc:= subs(ma=mc, norm=cpodom, eval(itera));

Décalage de la matrice initiale

> mc:= inverse( ma);

Convergence : la direction se stabilise

> (iterc@@30)(v0); iterc(%); vp3:= cpodom(evalm(mc &* %));

Augmentation de la précision

> Digits:=20: vm:= (iterc@@100)(v0); vm1:= iterc(%); vq1:= cpodom(evalm(mc &* %));
evalm(vm1-vm); norm(%,2); Digits:=10:

Vitesse : 11 décimales en 100 étapes

> (10.^11)^(1/100);

Vérification

> [eigenvals(map(evalf, mc))]; sort(map(abs,%)); %[-1]/%[-2];

La dernière valeurs propre de ma est donc

> mu4:= evalf(1/vq1, 20);

Et le vecteurs propres

> wp4:= pldx[xpurge](vm): Vector[column](wp4);
ma &* wp4 - 'mu4'*wp4; evalf(evalm(%), 20): norm(%,2);
# remarque sur les arrondis de calcul : on perd deux décimales sur le résultat

>

Calcul de la dernière valeur propre

On pourrait se contenter d'une soustraction...

Réécriture de l'itérateur

> iterd:= subs(ma=md, norm=cpodom, eval(itera));

Décalage de la matrice initiale

> ts:= -30: md:= inverse( ma-ts);

Convergence

> (iterd@@30)(v0); iterd(%); vp3:= cpodom(evalm(md &* %));

Augmentation de la précision

> Digits:=20: vm:= (iterd@@30)(v0); vm1:= iterd(%); vq2:= cpodom(evalm(md &* %));
evalm(vm1-vm); norm(%,2); Digits:=10:

Vitesse : 12 décimales en 30 étapes

> (10.^12)^(1/30);

Vérification

> [eigenvals(map(evalf, md))]; sort(map(abs,%)); %[-1]/%[-2];

La dernière valeurs propre de ma est donc

> mu3:= evalf(ts+1/vq2, 20);

Et le vecteurs propres

> wp3:= pldx[xpurge](vm): Vector[column](%);
ma &* wp3 - 'mu3'*wp3; evalf(evalm(%), 20): norm(%,2);
# remarque sur les arrondis de calcul : on perd deux décimales sur le résultat

>

Test général

> conj:= z-> map(conjugate,z);

> Digits:=20 : evalm(wp1), conj(wp1), evalm(wp3), evalm(wp4):
pas:= transpose(matrix([%])): 'pas' = evalf(%,8);

> inverse(pas): sap:=pldx[xpurge](%): 'sap'= evalf(%, 8);

Et la matrice se diagonalise

> evalm(sap &* ma &* pas): evalf(pldx[xpurge](%,9), 10); Digits:=10:

>

>