Sperimentazioni dell'algoritmo

Ecco alcuni esempi di utilizzo del metodo delle potenze:

Proviamo innanzitutto che l'algoritmo funziona: lo proviamo cioè su una matrice di cui conosciamo lo spettro, per esempio una matrice diagonale: gli autovalori sono infatti gli elementi sulla diagonale.

» A

A =

    20     0     0
     0     5     0
     0     0     3

» [l,i]=potenze(A, rand(3,1), 10e-15, 2000)

l =

  20.00000000000000


i =

    14
»

Chiaramente l'autovalore dominante era l'elemento $a_{1,1}$ con valore 10; l'algoritmo è riuscito a trovare questo elemento in 14 passi e con la precisione richiesta.

La prova successiva consiste in una matrice i cui elementi sono interi, ma scelti casualmente:

» A

A =

    45     2     9     0
     0     9     8    50
     1     2     3     4
     6     9     1     3

» [l,i]=potenze(A, rand(4,1), 10e-15, 2000)

l =

  46.05467285995372


i =

    55

» eig(A)

ans =

  46.05467285995431
 -15.19630313772274
  26.92844782249211
   2.21318245527628

»

Vista la casualità della matrice per avere una riprova che l'algoritmo trovasse effettivamente l'autovalore di modulo massimo, si è utilizzato il comando ``eig'' che Matlab fornisce per il calcolo dello spettro di una matrice; come si vede il metodo delle potenze trova una buona approssimazione dell'autovalore dominante. Si tenga anche presente che il numero di passi necessari per avere una soluzione e la precisione di questa dipendono dal vettore iniziale: essendo scelto casuale questi valori possono cambiare, a volte anche in modo rilevante.

Naturalmente il metodo delle potenze funziona anche in caso di matrici e vettori iniziali complessi, come si nota nel seguente esempio:

» A

A =

 10     0     0         0
  0     3     0         0
  0     0     0 + i     0
  0     0     0         0 - i

» [l,i]=potenze(A, rand(4,1)+i*rand(4,1), 10e-15, 2000)

l =

  10.00000000000000


i =

    16

»

Abbiamo provato un semplice esempio di una matrice diagonale con elementi complessi, ed anche il vettore casuale di partenza è stato scelto complesso.

Proviamo adesso un caso in cui l'algoritmo non converge:

» A

A =

     4   141  -576   432
     1     0     0     0
     0     1     0     0
     0     0     1     0

» [l,i]=potenze(A, rand(4,1), 10e-15, 2000)

l =

   3.09646771834061


i =

        2000

» eig(A)

ans =

 -12.00000000000000
  12.00000000000000
   3.00000000000000
   1.00000000000000

»

Osservando l'output del comando ``eig'' ci accorgiamo subito dove risiede il problema: ci sono due autovalori che hanno modulo massimo; in questo caso ``potenze'' impiega tutte le $2000$ iterazioni concesse senza raggiungere un'approssimazione accettabile. Questo è dovuto al fatto che durante la costruzione dell'algoritmo avevamo supposto che esistesse uno ed un solo autovalore dominante, qui ce ne sono due e l'algoritmo non riesce a convergere. Un'altra cosa da notare è la forma della matrice che abbiamo utilizzato per quest'ultimo test: una matrice siffatta si chiama matrice di compagnia e vedremo tra poco come si possono costruire matrici di questo tipo.