Matlab - Metodo di Jacobi

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

---

Torna alla pagina di Analisi Numerica