Calcul de valeurs propres

Pour une matrice A la valeur propre de plus grand module a un rôle dominant lorsqu'on calcule les puissances de A, de façon analogue au terme de plus haut degré d'un polynôme. A la limite A se comporte approximativement comme une multiplication scalaire par cette valeur propre. Cette remarque est utilisée pour calculer le vecteur propre correspondant ainsi que la valeur propre en itérant simplement la multiplication par A sur un vecteur donné.

1.Algorithme

Soit A une matrice $n\times n$ symétrique réelle. On note $\lambda_1,\lambda_2,\dots,\lambda_n$ ses valeurs propres (supposées réelles) et $v_1,v_2,\dots,v_n$ les vecteurs propres correspondants (normés).

On suppose que les $\lambda_i$ sont de modules distincts et on les ordonne par modules décroissants : $\vert\lambda_1\vert\gt\vert\lambda_2\vert\gt\dots\vert\gt\vert\lambda_n\vert$

Soit x⁽⁰⁾ un vecteur de Rⁿ et la suite x^(k) définie par la relation de récurrence

$$x^{(k+1)}={{\textstyle Ax^{(k)}}\over{\textstyle \Vert Ax^{(k)}\Vert}}$$

où $\Vert \ \Vert$ est la norme euclidienne sur Rⁿ.

Montrer que $x^{(k)}\rightarrow v_1$ quand $k\rightarrow\infty$.On forme alors $A^{(1)}=A-\lambda_1v_1v_1^t$(v₁^t vecteur-ligne transposé du vecteur-colonne v₁).

Montrer que les valeurs propres de A⁽¹⁾ sont $\lambda_2,\dots,\lambda_n,0$.

Que trouve-t-on si on forme une suite par la relation de récurrence :

$$x^{(k+1)}={{\textstyle A^{(1)}x^{(k)}}\over{\textstyle \Vert A^{(1)}x^{(k)}\Vert}}$$

On en déduit l'algorithme :

    pour i de 1 à n

        calcul de la limite v de la suite x^(k)

        calcul du $\lambda$ correspondant

        remplacement de A par A-$\lambda$vv^t

Comment proposez-vous de calculer $\lambda$ à partir du vecteur v ?

Quel critère peut permettre de vérifier numériquement que v est bien un vecteur propre ?

2.Programmation

On dispose d'un programme qui pour plusieurs matrices A calcule les valeurs propres et vecteurs propres d'après la méthode précédente.

On utilise les types :

Vecteur=ARRAY[1..MatMax] OF REAL

Matrice=ARRAY[1..MatMax,1..MatMax] OF REAL

Ecrire les procédures suivantes dont seul l'en-tête a été indiqué :

PROCEDURE CalculVecteur(A:Matrice;Dim:INTEGER;VAR V:Vecteur);

qui pour la matrice A de dimension Dim retourne le vecteur V limite de la suite x^(k) définie au 1. ; on prendra pour x⁽⁰⁾ un vecteur à composantes aléatoires (fonction Random) et on limitera k à la valeur ItMax, constante définie en début de programme.

FUNCTION ValeurPropre(A:Matrice;Dim:INTEGER;V:Vecteur):REAL

qui retourne la valeur propre correspondant au vecteur propre V de la matrice A

PROCEDURE Deflation(VAR A:Matrice;Dim:INTEGER;V:Vecteur;Val:REAL);

qui remplace la matrice A par $A-Val\times v\times v^t$(v vecteur de norme euclidienne 1).

On utilisera la fonction :

FUNCTION Norme(X:Vecteur;Dim:INTEGER):REAL;

déjà définie dans le programme.

Exécuter le programme pour les matrices proposées puis envisager une procédure supplémentaire pour vérifier les résultats.

$\bullet$ Partie théorique

On décompose x⁽⁰⁾ sur la base $v_1,v_2,\dots,v_n$ :

$$x^{(0)}=\sum_{i=1}^nc_iv_i$$

On considère aussi la suite y^(k) définie par

$$y^{(0)}=x^{(0)}\qquad y^{(k+1)}=Ay^{(k)}$$

On a

$$x^{(k+1)}={{\textstyle Ay^{(k)}}\over{\textstyle \Vert Ay^{(k)}\Vert}}$$

On a alors

$$Ay^{(0)}=\sum_{i=1}^nc_iAv_i =\sum_{i=1}^nc_i\lambda_iv_i$$

donc

$$y^{(1)}=\sum_{i=1}^nc_i\lambda_iv_i$$

et par récurrence, à condition que $c_1\ne 0$ :

$$y^{(k)}=\sum_{i=1}^nc_i\lambda_i^kv_i =c_1\lambda_1^k(v_1+ \... ..._1}}{{\textstyle \lambda_i^k}\over{\textstyle \lambda_1^k}}v_i)$$

$$\Vert y^{(k)}\Vert=\vert c_1\lambda_1^k\vert\left(1+ {{\text... ...e \lambda_i^{2k}}\over{\textstyle \lambda_1^{2k}}}\right)^{1/2}$$

Donc

$$x^{(k)}=\sum_{i=1}^n\alpha_i^kv_i$$

Comme $\vert\lambda_i\vert<\vert\lambda_1\vert$on a ${{\textstyle \lambda_i^k}\over{\textstyle \lambda_1^k}}\rightarrow 0$ et ainsi

$$\vert\alpha_1^k\vert\rightarrow 1 \qquad\alpha_i^k\rightarrow 0$$

Pour A⁽¹⁾ on a

$$A^{(1)}v_1=(A-\lambda_1v_1v_1^t)v_1=Av_1-\lambda_1v_1=0$$

$$A^{(1)}v_i=(A-\lambda_1v_1v_1^t)v_i=Av_i=\lambda_iv_i$$

donc A⁽¹⁾ a les mêmes vecteurs propres que A avec les valeurs propres $\lambda_2,\dots,\lambda_n,0$.

Si on forme alors une suite définie par la relation de récurrence :

$$x^{(k+1)}={{\textstyle A^{(1)}x^{(k)}}\over{\textstyle \Vert A^{(1)}x^{(k)}\Vert}}$$

on va trouver v, vecteur propre associé à la plus grande valeur propre de A⁽¹⁾ en module, c'est-à-dire $\lambda_2$.

Pour calculer $\lambda$ on peut simplement utiliser la définition $Av=\lambda v$ et prendre $\lambda={{\textstyle (Av)_k}\over{\textstyle v_k}}$ ; il faut s'assurer que $v_k\ne 0$ et que $\lambda$ est bien le même suivant k. Une vérification peut alors être faite en calculant $\Vert Av-\lambda v\Vert$.

Il est cependant préférable de calculer $\lambda={{\textstyle (Av\vert v)}\over{\textstyle (v\vert v)}}$, la colinéarité de v et Av étant vérifiée par le fait que ${{\textstyle \vert(Av\vert v)\vert}\over{\textstyle \Vert Av\Vert\Vert v\Vert}}$vaut 1.

$\bullet$ Programmation

FUNCTION Norme(X:Vecteur;Dim:INTEGER):REAL;
VAR I:INTEGER;
    S:REAL;
BEGIN
  S:=0.0;
  FOR I:=1 TO Dim DO S:=S+SQR(X[I]);
  Norme:=SQRT(S)
END;

PROCEDURE CalculVecteur(A:Matrice;Dim:INTEGER;VAR V:Vecteur);
VAR I,J,K:INTEGER;
    U,Up:Vecteur;
    N,T:REAL;
BEGIN
  FOR I:=1 TO Dim DO Up[I]:=random;
  K:=1;
  REPEAT
    U:=Up;
    FOR I:=1 TO Dim DO
    BEGIN
      T:=0;
      FOR J:=1 TO Dim DO T:=T+A[I,J]*U[J];
      Up[I]:=T
    END;
    N:=Norme(Up,Dim);
    IF Up[1]<0 THEN N:=-N;
    T:=0;
    FOR I:=1 TO Dim DO
    BEGIN
      Up[I]:=Up[I]/N;
      T:=T+ABS(Up[i]-U[I])
    END;
    K:=K+1;
  UNTIL (K>ItMax) OR (T<Eps);
  V:=Up
END;

FUNCTION ValeurPropre(A:Matrice;Dim:INTEGER;V:Vecteur):REAL;
VAR I,J:INTEGER;
    N,D:REAL;
BEGIN
  N:=0.0; D:=0.0;
  FOR I:=1 TO Dim DO
  BEGIN
    D:=D+V[I]*V[I];
    FOR J:=1 TO Dim DO N:=N+V[I]*A[I,J]*V[J];
  END;
  ValeurPropre:=N/D
END;

PROCEDURE Deflation(VAR A:Matrice;Dim:INTEGER;V:Vecteur;Val:REAL);

VAR I,J:INTEGER;

BEGIN
  FOR I:=1 TO Dim DO
  BEGIN
    FOR J:=1 TO Dim DO A[I,J]:=A[I,J]-val*V[I]*V[J];
  END;
END;

Pour des matrices de valeurs propres connues on constate que la qualité numérique des dernières valeurs propres est bien moins bonne que celle des premières. Ceci est du à l'accumulation d'erreurs au cours des calculs successifs de A⁽¹⁾,A⁽²⁾,...,A^(n-1). En pratique cette méthode, de mise en oeuvre assez simple, ne peut être valablement utilisée que pour calculer les quelques valeurs propres de plus grand module d'une matrice. Ceci est justement recherché en mécanique (modes fondamentaux de vibration) ou en chimie (niveaux d'énergie).