|
|
|
|
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
|
|
|
I should add that I am using a zero starting guess, and the convergance test is testing the relative residual.
|
|
|
On Sun, Oct 28, 2012 at 5:40 PM, Joza < [hidden email]> wrote:
Your 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
|
|
|
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!
|
|
|
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!!??
|
|
|
----- 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
|
|
|
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
|
|
|
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
|
|
|
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
|
|
|
Since I permute A's rows using dmperm, how do I also permute b in the same way?
|
|
|
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
|
|
|
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!!!
|
|
|
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)
|
|
|
On Tue, Oct 30, 2012 at 6:10 AM, Juan Pablo Carbajal
< [hidden email]> wrote:
Juan,
Great pointer to agora.octave.org. I hadn't heard of it before.
I just posted an update to your script that attempts to create a
better permutation matrix. It seems to me that to make a matrix
diagonal dominant, it is a necessary condition that the largest
magnitude value in any row be on the diagonal. I added an initial
permutation "guess" p that is simply the indices of max(abs(A)). I
then test to see if all of the rows are unique, and, if not fall back
to using dmperm.
This will definitely help with the test case that the initial poster was using.
Bill
_______________________________________________
Help-octave mailing list
[hidden email]
https://mailman.cae.wisc.edu/listinfo/help-octave
|
|
Locked
