Re: Implementing a Jacobi iterative method for Ax=b

[Joza]

So I've taken all this on board, and finally got an algorithm which as far as
I am concerned SHOULD work...

Please note that I'm using the "correction form" of the Jacobi algorithm:
x_new = x_old + D_inv*(b - A*x_old);
where A is from Ax=b (but possibly permuted).

You can derive it from the usual x_new = D_inv*( b - L*x_old - U*x_old )
where A = D + L + U.

So here is the code:
************************************************************************************************************************************************
function [x k] = Jacobi(A, b, tol)

k = 0;
n = size(A,1);
x_old = zeros(n,1);
converged = 0;

W = A(dmperm(A), :); % dmperm permutes rows of A so that no diagonal
elements are zero.
v = b(dmperm(A),1);  % must permute b the same way so that x is unchanged.

for i=1:n
        D_inv(i,i) = 1/W(i,i);
end

while ~converged

                x_new = x_old + D_inv*(v - W*x_old);
                k = k + 1;
                
                if norm(v - W*x_new, inf)/norm(v, inf) < tol
                        x = x_new;
                        converged = 1;
                else
                        x_old = x_new;
                end
        
end
*******************************************************************************************************************8
I have been given a starting A and b and told to use the zero starting
approximation.

When I run the algorithm, I do not get a return x or k, all the terminal
gives me is ans = [NaN, NaN, NaN ...]

So clearly this is running into difficulty somewhere and not finishing
correctly, but I really can't see where. Everything should work as far as I
can see!!!



--
View this message in context: 
http://octave.1599824.n4.nabble.com/Implementing-a-Jacobi-iterative-method-for-Ax-b-tp4645833p4645866.html
Sent from the Octave - General mailing list archive at Nabble.com.