Torna alla pagina di Analisi Numerica
:: Matlab - Metodo di Jacobi ::
% Metodo iterativo di Jacobi
% --------------------------
% restituisce la soluzione del sistema x ed il numero
% di iterazioni effettuate N
% 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]=jacobi(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 le variabili xold e x che utilizzerò e
% aggiornerò nelle varie iterazioni
xold = x0;
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 vecchi valori
% della variabile x
for i=1:n
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) * xold(j);
end
for j=i+1:n,
s = s + a(i,j) * xold(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
res = norm(r);
% termino l'iterazione assegnando alla variabile xold
% la x attuale, così che il metodo possa proseguire
xold = x;
end