Ricerca autovalori con metodi delle potenze

[LLLorenzzz]

Ciao a tutti

sto svolgendo un'esercitazione per un corso di analisi numerica con matlab.
Mi si chiede di implementare funzioni che calcolino gli autovalori di una matrice data con i metodi delle potenze in norma uniforme, euclidea e col metodo di Wielandt (potenze inverse con shift, che ho implementato in norma euclidea).
Mentre i metodi "diretti" cioè quelli che cercano l'autovalore massimo convergono benissimo, il metodo di Wielandt, che cerca l'autovalore vicino ad un'approssimazione data, è lentissimo e non converge praticamente mai (usando anche 1000 iterazioni e una tolleranza di $10^-6$). Per guadagnare velocità ho provato anche a fattorizzare (LU) la matrice prima di iniziare il metodo, così da abbassare il costo a $n^2$ anziché $n^3$(si tratterebbe di risolvere un sistema lineare a ogni passo altrimenti). Volevo chiedervi se secondo voi questo è normale e da cosa può dipendere.

grazie mille ciaaaao

[alvinlee88]

A me non sembra normale, credo possa dipendere da quanto è distante dal più vicino autovalore la tua approssimazione. In generale un'approssimazione \( \displaystyle \mu_k \) di \( \displaystyle \lambda_k \) autovalore è "buona" se \( \displaystyle min_{i=1,...n}|\lambda_i-\mu_k|=|\lambda_k-\mu_k|\ne0 \) .
Se questa condizione è rispettata, esiste l'autovalore dominante di \( \displaystyle (A-\mu_kI)^{-1} \) , ossia \( \displaystyle \frac{1}{\lambda_k-\mu_k} \) , quindi applicando il metodo diretto a \( \displaystyle (A-\mu_kI)^{-1} \) esso dovrebbe convergere (a quanto ne so in soli 3-4 passi, anche se non l'ho mai sperimentato numericamente) all'autovettore di \( \displaystyle (A-\mu_kI)^{-1} \) relativo a \( \displaystyle \frac{1}{\lambda_k-\mu_k} \) ,che coincide con l'autovettore di \( \displaystyle A \) relativo a \( \displaystyle \lambda_k \) .
Il costo per trovare questo autovettore (fattorizzando inizialmente la matrice) è \( \displaystyle O(n) \) per matrici simmetriche, \( \displaystyle O(n^2) \) altrimenti.

Ciao

[LLLorenzzz]

Grazie mille per aver risposto!

ci sto lavorando e in effetti ho ottenuto dei miglioramenti...può dipendere dalla norma che uso? se lo implemento in norma uniforme non fa mai più di una ventina di passi anche su matrici belle grandi.

grazie ancora ciao!

[alvinlee88]

Nel corso di calcolo scientifico che ho seguito, in cui mi è stato spiegato questo metodo, veniva usata appunto la norma infinito, o uniforme, e i discorsi sulla convergenza fatti sopra erano appunto usando questa norma.

Se non ci hai già pensato, potresti prima di far partire il metodo lridurre a tua matrice $A$ in forma di hessemberg superiore (o forma tridiagonale se è simmmetrica), e comunque guarda se riesci a verificare se le tue approssimazioni sono "buone" nel senso definito sopra. Se lo sono a quanto so dovrebbe metterci un pò meno di 20 passi...fammi sapere, ciao.

[LLLorenzzz]

Ciao alvinlee88!
ho fatto come mi hai consigliato ed effettivamente è molto più veloce adesso
ti ringrazio davvero ciao!

[alvinlee88]

Di niente, figurati
Comunque il ridurre preventivamente le matrici in forma di hessemberg superiore o tridiagonale e' un fatto usato spesso nella pratica, poichè in genere uno stesso algoritmo funziona molto più velocemente su una matrice di questa forma.

P.S. A mio parere un ottimo libro per l'analisi numerica (e più nello specifico l'algebra lineare numerica) è "Applied Numerical Linear Algebra", di J.W. Demmel, molto completo sia a livello della matematica richiesta dal settore, sia per quanto riguarda la pratica (software, sperimentazione numerica, dimostrazioni anche molto tecniche riguardo stabilità e condizionamento dei problemi ecc.)
L'ho usato poco (come integrazione a un corso che ho seguito) ma mi è sembrato proprio fatto bene.
ciao

[LLLorenzzz]

bene allora vedrò di procurarmelo visto che il mio professore non ha dato nessun testo per il corso!
Grazie mille di tutto alvinlee! grazie davvero!

[MathematicianGirl]

Uffi! Io non riesco a far girare il programma...ci sto impazzendo! Qualcuno può aiutarmi con l'implementazione del metodo delle potenze?
Riesco sempre a capire bene gli algoritmi, ma quando si tratta di implementarli vado in tilt.....
Vi ringrazio in anticipo!
P.S. Qua di seguito vedete il mio programmino....magari sareste così gentili da aiutarmi a capire cosa c'è che non va.
GRAZIE 1000!

[m,n]=size(A);
if m~=n
disp('La matrice deve essere quadrata')
return
end
k=0;
zold=y0;
told=zold/norm(zold,inf);
while k<nmax
znew=A*told;
tnew=znew/norm(znew,inf);
told=tnew;
snew=(ctranspose(tnew)*znew)/(ctranspose(tnew)*tnew);
k=k+1;
end
snew;