[Math] Jacobi Iteration Algorithm in MATLAB

Question

So I have to write a Matlab algorithm to perform a Jacobi iteration. It needs to be executed as >jacobi(A, b, x0, tol, Niter). Here is my algorithm so far:

function x1 = myjacobi(A, b, x0, tol, Niter)
%Step 1
k = 1;
n = length(b);
%Step 2
while(k < Niter)
    %Step 3
    for j = 1:n
        x(j) = ((b(j) - A(j,[1:j-1,j+1:n]) * x0([1:j-1,j+1:n])) / A(j,j));
    end
    %Step 4
    if norm(x0) < tol
        for i=1:n
            disp(x(i));
        end
    end
    %Step 5
    k = k+1;
    %Step 6
    for i=1:n
        x0(i) = x(i);
    end   
end
%Step 7
disp ('Maximum number of iterations exceeded');
end

However my program should display a message that Jacobi failed to converge if the method doesn't converge after Niter maximum iterations. And it needs to compute a solution of Ax=b starting at x0=0 and stop when either norm(x(k) – x(k-1))/(x(k)) < tol or k = Niter. It also needs to compute the residual r = norm(Ax(k) – b) and the error e = norm(x(k) – x*). I'm not sure where to put these components in my code. Any assistance would be appreciated.

Answer 1

I think this is what you need:

function x = myjacobi(A, b, x0, tol, Niter)

% Step 1
k = 1;
n = length(b);
x = x0;
D = diag(diag(A));
R = A - D;
Dinv = inv(D);

% Step 2
while (k < Niter)
    % Step 3
    xold = x;
    x = Dinv * (b - R*xold);

    % Step 4
    if norm(x - xold) / norm(x) < tol
        disp('Stop iterating because norm(x(k)-x(k-1))/norm(x(k)) < tol')
        break;
    end

    % Step 5
    k = k+1;
end

% Step 6
if k == Niter
    disp('Maximum number of iterations exceeded');
end

% Compute residual
r = norm(A * x - b);
fprintf('Residual: %.3e\n', r);

% Compute error
xstar = A\b;
e = norm(x - xstar);
fprintf('Error: %.3e\n', e);