2
$\begingroup$

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.

$\endgroup$
1
$\begingroup$

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);
$\endgroup$
  • $\begingroup$ Thank you very much, that improved things! Could you explain your logic for the while loop? Could I have done the same thing with my while and for loop? $\endgroup$ – Harshi Nov 2 '17 at 6:58
  • $\begingroup$ Actually, I used your while loop as starting point, so yes you could have done the same thing with your while loop :). I replaced your 'step 3' with a shorter notation, by making use of matrix multiplications. Furthermore, I think you checked for the wrong thing at 'step 4', so I changed that too. What you do at 'step 6': there is no need to do this in a for loop. Just saying 'x0=x' would do the same. I do this at 'step 3'. $\endgroup$ – EdG Nov 2 '17 at 7:10
  • $\begingroup$ That makes sense, thank you! But n is never used, and it needs to change as n gets larger. $\endgroup$ – Harshi Nov 2 '17 at 20:07
  • $\begingroup$ You are right in the sense that n is never used. So that specific line could be removed. I don't understand, however, why we should use it as n gets larger. I don't see n appearing in your formulas, so could you explan that? $\endgroup$ – EdG Nov 2 '17 at 23:55

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.