Come plottare in matlab

[Studente Anonimo]

Mi viene richiesto di implementare il Metodo del gradiente con precondizionamento e di plottare la norma di \( r_k \) vs numero di iterazioni \(k\). Ora. Io ho una funzione my_gradient al cui interno ho un ciclo while, come faccio nella main a tirare fuori l'informazione della norma di \(r_k\) ad ogni passo alla \(k\)-esima iterazione e poi a plottare il grafico richiesto?

[feddy]

Ad ogni iterazione $k$ calcoli il $||r_k||$ e lo salvi in un vettore $r$. E' sufficiente che la tua funzione my_gradient ritorni il vettore $r$, e poi puoi semplicemente plottare in scala semilogaritmica l'andamento della norma. Se non ti è chiaro posta direttamente il codice.

Mi pare che su matlab la sintassi per i valori di ritorno di una funzione siano tipo

[val1, val2] = function(par1,par2,...)

[Studente Anonimo]

Questo è il mio codice

function [out] = my_gradient(A,P,b,x0,max,tol)
xk=x0;
k =1 ;
while ( ( norm(b - A * xk)./norm(b)  >= tol ) && k <= max )
     rk = b - A * xk;
     yk = P \ rk;
     alphak = yk' * (rk) / ((A * yk)' * yk);
     xk = xk + alphak * yk;
    
end
if k == max
    warning('Does not convergence in max iteration');
end
out = xk;
end

Il punto è che restituisco la "soluzione" del sistema lineare \( Ax = b \), e quindi non so dare un output in più e poi come dire nella main che mi interessa solo uno degli output.

[feddy]

Ti basta una cosa del tipo

function [out,rk_hist] = my_gradient(A,P,b,x0,max,tol)
 

dove

rk_hist

altro non è che un vettore che ad ogni passo k si salva rk
così la funzione my_gradient ritorna due cose: la prima è l'array che contiene la soluzione una volta che il metodo è arrivato a convergenza. La seconda è un vettore che contiene la storia dei residui che nel while cambia dimensione ad ogni iterata

rk_hist = [rk_hist, b-A*xk];

dovrebbe andare. Per dire nel main che ti interessa solamente un output e il resto no, basta scrivere

[out,~] = my_gradient(A,P,b,x0,max,tol)

(un po' simile a Python se lo conosci)

ma comunque in realtà dovresti proprio scrivere nel main

[out,rk_hist] = my_gradient(A,P,b,x0,max,tol);

siccome tu non vuoi solo la soluzione, ma anche la storia dei residui. Poi ti basta fare

semilogy(1:numel(rk_hist),rk_hist,'*');

mi pare che il comando per restituire il numero di elementi del vettore sia

numel

, sto andando a memoria, non uso MatLab da tanto.

[feddy]

Tra l'altro ho notato che nel tuo codice ci sono un po' di errori: non incrementi k, e poi l'algoritmo mi sembra sia scritto male... da dove stai studiando per implementarlo?

Su wiki c'è una possibile versione in pseudocodice: https://en.wikipedia.org/wiki/Conjugate ... ent_method

[Studente Anonimo]

Grazie mille!
Mah... a dire il vero lo sto studiando "a caso" nel senso seguendo il corso di analisi numerica il prof ci ha fatto una lezione introduttive di 2 ore sui comandi base di matlab e poi basta, si aspetta che lo sappiamo usare... -.-È vero che si concentra più sull'aspetto teorico: dimostrare perché un metodo piuttosto che un altro funziona, ma tutti gli esercizi dopo l'aspetto teorico contengono una parte di applicazione. Io cerco di volta in volta come fare una cosa quando mi serve...
Non conosco python ma un po' di c++, anche se è due anni che non programmo più in c++

[feddy]

C++ FTW

Se conosci C++, Python (e soprattutto MatLab) non ti daranno problemi in termini di difficoltà di apprendimento.

Tornando al problema, se vuoi provare ad implementarlo basandoti sul listato di Wiki poi posso darti un check. Non credo tu possa avere particolari difficoltà.

[feddy]

Intendo implementarlo in MatLab. In C++ dovresti usare librerie esterne come Eigen ed è inutilmente complicato in questo caso

[Studente Anonimo]

Sono con il cellulare in metro. Quindi sorry se ho fatto qualche typo. Solo per la concezione del codice senza star lì a pensare a salvare la norma di rk fare così:

function [out] = my_gradient(A,P,b,x0,max,tol)
xk=x0;
rk=b-A*xk;
zk=P\rk;

k =1 ;
while ( ( norm(b - A * xk)./norm(b)  >= tol ) && k <= max )
     alphak=dot(rk' * zk)/dot(zk' *A*zk);
     xk= xk+ alphak*zk;
     rk=rk-alphak*A*zk;
     k=k+1;
    
end
if k == max
    warning('Does not convergence in max iteration');
end
out = xk;
end

Ps la condizione di uscire dal while la vuole il prof così

[feddy]

Manca l'update di $z_k$, $p_k$,per il quale ti serve $\beta_k$. (sto facendo riferimento alla pagina wiki). Fai pure con calma

[feddy]

Ho trovato MatLab sul computer dell'ufficio, quindi ho scritto quello che mancava del tuo codice. L'ho testato al volo e funziona sulla matrice che discretizza $-\Delta u = 1$ sul quadrato. E' una matrice notoriamente simmetrica e definita positiva. Sotto l'andamento della norma del residuo.

clear all
close all

N = 15;
A = gallery("poisson",N); %matrice N*N, corrispondente alla discr. del problema di Poisson in 2D
b = ones(N^2,1);
P = eye(N^2);
x0 = ones(N^2,1);
max = 200;
tol = 1e-12;
[out,rk_hist] = my_gradient(A,P,b,x0,max,tol);
semilogy(1:numel(rk_hist),rk_hist,'*')
xlabel('iteration PCG')
ylabel('Residual norm')

function [out,rk_hist] = my_gradient(A,P,b,x0,max,tol)
xk=x0;
k=0;
rk = b-A*xk;
zk = P\rk;
pk = zk;
rk_hist = [norm(rk)];
while ( ( norm(rk)/norm(b)  >= tol ) && k <= max )
    Apk = A*pk;
    alphak = (rk'*zk)/(pk'*(Apk));
    xk = xk + alphak*pk;
    rkp1 = rk- alphak*(Apk);
    rk_hist = [rk_hist,norm(rkp1)];
    if (norm(rkp1,2)<tol)
        out = xk;
        break;
    else
        zkp1 = P\(rkp1);
        betak = (rkp1'*zkp1)/(rk'*zk);
        pk = zkp1 + betak*pk;
        rk=rkp1;
        zk = zkp1;
        k = k+1;
    end
    
end

if k == max
    warning('Does not convergence in max iteration');
end

out = xk;
end

[Studente Anonimo]

Scusa la risposta tardiva,
Grazie mille!

[feddy]

Di nulla