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.
- Codice:
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
