# Implementing a Jacobi iterative method for Ax=b Classic List Threaded 15 messages Open this post in threaded view
|

## Implementing a Jacobi iterative method for Ax=b

 Hi guys. I am trying to implement an iterative method for solving Ax=b using the Jacobi method. When I run the code with any tolerance it doesn't end, however. Any ideas why? The matrix A is relatively sparse. Should this mean I should do row permutations on it before the iteration begins? ********************************************************* function [x iter] = Jacobi(A, b, tolerance) iter = 0; n = size(A,1); x_old = zeros(n,1); converged = 0; M = diag(diag(A)); N = M - A; while ~converged                 x_new = M\(N*x_old + b);                 iter = iter + 1;                                 if norm(b - A*x_new)/norm(b) < tolerance                         x = x_new;                         converged = 1;                 else                         x_old = x_new;                 end         end
Open this post in threaded view
|

## Re: Implementing a Jacobi iterative method for Ax=b

 I should add that I am using a zero starting guess, and the convergance test is testing the relative residual.
Open this post in threaded view
|

## Re: Implementing a Jacobi iterative method for Ax=b

 On Sun, Oct 28, 2012 at 5:40 PM, Joza <[hidden email]> wrote: > I should add that I am using a zero starting guess, and the convergance test > is testing the relative residual. > > > > -- > View this message in context: http://octave.1599824.n4.nabble.com/Implementing-a-Jacobi-iterative-method-for-Ax-b-tp4645833p4645834.html> Sent from the Octave - General mailing list archive at Nabble.com. > _______________________________________________ > Help-octave mailing list > [hidden email] > https://mailman.cae.wisc.edu/listinfo/help-octaveYour code has one error and one unnecessary calculation. Error: The step of the Jacobi iteration is x_new = M\(b - N*x_old); Unncesessary calculations M is a diagonal matrix, therefore the inverse is just the inverse of each diagonal element (take care of zero elements!). Improve the speed of you code as follows M = diag(A); do not build the matrix, just keep the vector x_new = (b - N*x_old)./M; The method is not always convergent. Check that the input matrix guarantees convergence for your tests. For real life usage, limit the maximum numbe rof interations. When aborting because of max number of iteration provde a warning and return enough information for the user to take a decision. http://en.wikipedia.org/wiki/Jacobi_method#ConvergenceHope that helps _______________________________________________ Help-octave mailing list [hidden email] https://mailman.cae.wisc.edu/listinfo/help-octave
Open this post in threaded view
|

## Re: Implementing a Jacobi iterative method for Ax=b

 Thanks for the input! My method was correct, but you have informed me that I should be using the "correction form". How can I permute an arbitrary matrix so that it's diagonal entries are not zero? this will affect my solution I'm sure!
Open this post in threaded view
|

## Re: Implementing a Jacobi iterative method for Ax=b

 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!!??
Open this post in threaded view
|

## Re: Implementing a Jacobi iterative method for Ax=b

 ----- Original Message ----- > From: Joza <[hidden email]> > To: [hidden email] > Cc: > Sent: Sunday, October 28, 2012 11:48 PM > Subject: Re: Implementing a Jacobi iterative method for Ax=b > > 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!!?? > "Whats going on!!??" - you forgot to put a semicolon at the end of         W = A(dmperm(A), :)             x = x_new     D_inv(i,i) = 1/W(i,i) lines. Regards,   Sergei. _______________________________________________ Help-octave mailing list [hidden email] https://mailman.cae.wisc.edu/listinfo/help-octave
Open this post in threaded view
|

## Re: Implementing a Jacobi iterative method for Ax=b

 In reply to this post by Joza On Sun, Oct 28, 2012 at 4:48 PM, Joza <[hidden email]> 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 _______________________________________________ Help-octave mailing list [hidden email] https://mailman.cae.wisc.edu/listinfo/help-octave
Open this post in threaded view
|

## Re: Implementing a Jacobi iterative method for Ax=b

 In reply to this post by Juan Pablo Carbajal-2 On 28 Oct 2012, at 17:59, Juan Pablo Carbajal wrote: > do not build the matrix, just keep the vector > x_new = (b - N*x_old)./M; Octave is clever enough to do that itself, diagonal matrices are a special class: >> typeinfo (diag (diag (M))) ans = diagonal matrix and have their own methods. c. _______________________________________________ Help-octave mailing list [hidden email] https://mailman.cae.wisc.edu/listinfo/help-octave
Open this post in threaded view
|

## Re: Implementing a Jacobi iterative method for Ax=b

 Oh, and if you permute A, I think you will also need to permute b so that they match, and then permute the result back to the original form.  I didn't try a matrix with any 0s on the diagonal. I used: A=[-500, 2, 3; 3, 200, 1; 4,5,600] b=[1;2;3] _______________________________________________ Help-octave mailing list [hidden email] https://mailman.cae.wisc.edu/listinfo/help-octave
Open this post in threaded view
|

## Re: Implementing a Jacobi iterative method for Ax=b

 Since I permute A's rows using dmperm, how do I also permute b in the same way?
Open this post in threaded view
|

## Re: Implementing a Jacobi iterative method for Ax=b

 On Sun, Oct 28, 2012 at 6:58 PM, Joza <[hidden email]> wrote: > Since I permute A's rows using dmperm, how do I also permute b in the same > way? > > > > -- > View this message in context: http://octave.1599824.n4.nabble.com/Implementing-a-Jacobi-iterative-method-for-Ax-b-tp4645833p4645857.html> Sent from the Octave - General mailing list archive at Nabble.com. It has been a long time since I did anything like this, but since you are only permuting the rows of A, I think you only need to permute the rows of b in the same way, and x will still be ok.  If you were permuting both rows and columns of A, you would need to permute the rows of b in the same way that you permute the rows of A, and "unpermute" the x vector to undo the permutation of the columns of A. That is only from memory and a little Sunday evening thinking, so I could have that mixed up. Since you only need to permute b, you could, where you have:  W = A(dmperm(A), :); You could put: P=dmperm(A); W=A(P,:); bp=b(p); And then do your calculations with bp instead of b _______________________________________________ Help-octave mailing list [hidden email] https://mailman.cae.wisc.edu/listinfo/help-octave
Open this post in threaded view
|

## Re: Implementing a Jacobi iterative method for Ax=b

 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!!!
Open this post in threaded view
|

## Re: Implementing a Jacobi iterative method for Ax=b

 Here are A and b if you want to test it out! A = [ 3 12 0 -1 0 0; 4 0 31 1 0 0; 2 1 0 0 17 -3; 27 2 0 0 0 1; 0 0 0 -1 1 11; 0 0 0 24 -1 0 ] (6x6) b = [ 39 117 12 98 14 55]^T (6x1)