Torna alla pagina di Analisi Numerica
:: Matlab - Metodo di Gauss-Seidel ::
% Metodo iterativo di Gauss-Seidel
% --------------------------------
% restituisce la soluzione del sistema x ed il numero
% di iterazioni effettuate it
% richiede come parametri di ingresso la matrice del sistema a,
% il termine noto b, il vettore iniziale x0, l'errore tollerato
% epsilon, il numero massimo di iterazioni consentite
function [x,it]=gaussseidel(a,b,x0,epsilon,nmax);
% calcola la dimensione della matrice a
[n,n] = size(a);
% inizializzo la variabile it per il conteggio delle iterazioni
it = 0;
% per verificare che l'errore sia inferiore a epsilon dobbiamo
% stimarlo a posteriori, a partire cioè da una quantità già
% calcolata. Uno stimatore è il residuo, definito da:
r = b - a * x0;
% posso utilizzare la norma del residuo per effettuare
% successivamente il test di arresto
res = norm(r);
% inizializzo la variabile x che utilizzerò e aggiornerò nelle
% varie iterazioni
x = x0;
% entro nel ciclo in cui avviene il calcolo delle soluzioni
% esco dal ciclo solo quando la norma del residuo è minore
% della tolleranza epsilon (1), o quando oltrepasso il
% numero massimo di iterazioni consentite (2)
% (1) (2)
% ------------- ---------
while res > epsilon & it < nmax
% incremento di 1 il numero di iterazioni
it = it + 1;
% ciclo in cui costruisco la sommatoria dei prodotti
% tra gli elementi della matrice e i valori già
% calcolati della variabile x
for i=1:n
% inizializzo la variabile s in cui andrò a formare
% la sommatoria
s = 0;
% devo spezzare il ciclo for in due distinti perché
% nella sommatoria non deve figurare il caso in
% cui j = i . Il primo for infatti arriva fino a
% (i-1), mentre il secondo parte da (i+1)
for j=1:i-1,
s = s + a(i,j) * x(j);
end
for j=i+1:n,
s = s + a(i,j) * x(j);
end
% costruzione dell'equazione al passo i-esimo
x(i) = (b(i) - s) / a(i,i);
end
% ricalcolo il residuo ad ogni iterata...
r = b - a * x;
% ...e subito dopo la sua norma, così da attualizzare
% il test di arresto che costituisce la condizione
% del ciclo
res = norm(r);
end