Re: Implementing a Jacobi iterative method for Ax=b

[Ozzy Lash]

On Sun, Oct 28, 2012 at 4:48 PM, Joza <address@hidden> wrote:
> Here is my latest code:
>
> function [x, k] = Jacobi2(A, b, tol)
>
> k = 0;
> n = size(A,1);
> x_old = zeros(n,1);
> converged = 0;
>
> for i=1:n
>         if A(i,i) == 0.0
>                 W = A(dmperm(A), :)
>                 break
>         else W = A;
>         end
> end
>
> for i=1:n
>         D_inv(i,i) = 1/W(i,i)
> end
>
> while ~converged
>
>                 x_new = x_old + D_inv*(b - W*x_old)   % Correction form
>                 k = k + 1;
>
>                 if norm(b - W*x_new)/norm(b) < tol
>                         x = x_new
>                         converged = 1;
>                 else
>                         x_old = x_new;
>                 end
> end
> ************************************************
> I reorder the matrix using dmperm so that no diagonals are zero.
>
> When i run this get the error:
>
> warning: broken pipe -- some output may be lost
>
> Whats going on!!??
>
>

I'm guessing that the message comes because you had to hit control-c to stop it.

A few things:

The way you have left it, the W matrix still has the diagonal
elements.  Also your update is not quite rite, it isn't

x_new = x_old + D_inv*(b - W*x_old)

it should be:

x_new = D_inv*(b - W*x_old)

assuming you remove the diagonal elements of W.  I changed your script
to the following, and it seems to work for a matrix I tried:

function [x, k] = Jacobi2(A, b, tol)

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

for i=1:n
        if A(i,i) == 0.0
                W = A(dmperm(A), :)
                break
        else W = A;
        end
end
R=W;
for i=1:n
        D_inv(i,i) = 1/W(i,i)
        W(i,i)=0
end

while ~converged

                x_new = D_inv*(b - W*x_old)   % Correction form
                k = k + 1;


                if norm(b - R*x_new)/norm(b) < tol
                        x = x_new
                        converged = 1;
                else
                        x_old = x_new;
                end
end