i was calculate using gauss jacobi method with 4x4 matrice, why does NaN appear? what should i do?
Show older comments
function[Xi,g,H]=jacobi(A,B,X,N)
clc;
A=[1 3 1 -1; 2 0 1 1; 0 -1 4 1; 0 1 1 -5]
B=[3;1;6;16];
X=[0;0;0;0];
Xi=inv(A)*B;
N=10;
abcd=X';
n=4;
X1=X;
for k=1:N,
for i=1:n,
S=B(i)-A(i,[1:i-1,i+1:n])*X([1:i-1,i+1:n]);
X1(i)=S/A(i,i);
end
g=abs(X1(i)-X(i))/X1(i)*100;
eror=norm(g);
X=X1;
abcd=[abcd;X'];
end
disp('abcd');
disp([abcd]);
disp('eror');
disp([eror]);
Answers (2)
A(2,2) equals 0, but you divide by 0.
What is H ?
clc;
A=[1 3 1 -1; 2 0 1 1; 0 -1 4 1; 0 1 1 -5];
A=[A(2,:);A(1,:);A(3,:);A(4,:)];
B=[3;1;6;16];
B=[B(2);B(1);B(3);B(4)];
X=[0;0;0;0];
N=20;
[Xi,g,H]=jacobi(A,B,X,N)
function[Xi,g,H]=jacobi(A,B,X,N)
Xi=inv(A)*B;
abcd=X';
n=4;
X1=X;
for k=1:N,
for i=1:n,
S=B(i)-A(i,[1:i-1,i+1:n])*X([1:i-1,i+1:n]);
X1(i)=S/A(i,i);
end
g=abs(X1(i)-X(i))/X1(i)*100;
eror=norm(g);
X=X1;
abcd=[abcd;X'];
end
disp('abcd');
disp([abcd]);
disp('eror');
disp([eror]);
end
John D'Errico
on 29 Oct 2023
Edited: John D'Errico
on 29 Oct 2023
You write a function that allows you to pass in A and B, but then you define A anb B in the function? Sigh. I think you do not understand why functions exist. Anyway...
A=[1 3 1 -1; 2 0 1 1; 0 -1 4 1; 0 1 1 -5]
Next, you want to read about the Jacobi method.
Somewhere along the way, your teacher should have said the Jacobi method ONLY applies to strictly diagonally dominant matrices. We can find this statement from that Wiki link I have copied here.
"A sufficient (but not necessary) condition for the method to converge is that the matrix A is strictly or irreducibly diagonally dominant. Strict row diagonal dominance means that for each row, the absolute value of the diagonal term is greater than the sum of absolute values of other terms:"
What does that mean? A matrix with any zero element on the diagonal can NEVER be convergent using Jacobi iteration. That is because you will be always dividing by zero. Those Nans are evidence of what happens.
So what should you do? Get a different matrix. Or, get a different algorithm, since Jacobi cannot be used to solve this problem.
5 Comments
Torsten
on 29 Oct 2023
After row permutation, the Jacobi method can be applied.
John D'Errico
on 29 Oct 2023
Edited: John D'Errico
on 29 Oct 2023
You are flat out incorrect in this.
Even if a row permutation would suffice on SOME matrices, It would no longer be the Jacobi method. Regardless, row permutations may not always be sufficient for Jacobi iterations on some matrices to be diagonally dominant, and therefore be convergent. In fact, it is easy enough to create a simple matrix where no row permutation is sufficient to restore diagonal dominance.
A = reshape(1:9,3,3)'
That matrix happens to be singular of course, but we can trivially change it to be not singular., yet still not be diagonally dominant.
rank(A)
A(3,3) = 9 - 100*eps(9);
rank(A)
Even a broader change will suffice for failure of the Jacobi method.
A(3,3) = 8;
rank(A)
I've written a simple Jacobi solver below, as jaciter. First, see that it works for a diagonally dominant problem.
A2 = A + 10*eye(3)
b = [2;3;5];
A2\b % known solution
[X,iter,del] = jaciter(A2,b,1e-14,1000)
So it converged to the known solution, after 76 iterations. Yet on problem A, it diverges
A\b % known solution
[X,iter,del] = jaciter(A,b,1e-14,1000)
See that we got a NaN here, probably from a difference like inf-inf.
Again though, on the matrix A, there is no set of row permutations that will restore diagonal dominance, even though A is not singular.
So, had you claimed that row permutations may, under some circumstances, improve convergence for a Jacobi iteration, you would be correct to some extent. There are SOME matrices where row permutation would help.
For example, change the problem by a simple set of row permutations...
A3 = A2([2 3 1],:)
b3 = b([2 3 1])
This does not change the system.
A3\b3
Yet now jaciter should fail.
[X,iter,del] = jaciter(A3,b3,1e-14,1000)
So it is true that on the problem A3\b3, row permutations will suffice to improve the problem, now allowing convergence.
function [X,iter,del] = jaciter(A,b,tol,itmax)
n = size(A,1);
D = diag(A);
LU = A - diag(D);
Dinv = diag(1./D);
X0 = rand(n,1);
del = inf;
iter = 0;
while (del>tol) && (iter <= itmax)
iter = iter + 1;
if iter >= itmax
warning('Maximum iterations exceeded')
end
X = Dinv*(b - LU*X0);
del = norm(X - X0);
X0 = X;
end
end
Torsten
on 29 Oct 2023
After row permutation, the Jacobi method can be applied.
Well, I meant in the actual case, of course.
Even on the actual problem, it is not strictly diagonally dominant, but yes, it is sufficient, on that specific problem to use a row permutation.
The problem is when you suggest that row permutations are sufficient to solve a problem, but NOT explain enough about the issues that you will now cause further issues down the line. Others might actually believe you, and think you have made a valid general claim because you have a good reputation on the site.
A=[1 3 1 -1; 2 0 1 1; 0 -1 4 1; 0 1 1 -5]
B=[3;1;6;16];
P = [2 1 3 4]; % a row permutation. It must be applied to both A and B.
A(P,:)
A you can see, now the diagonal elements are >= the sum of absolute values of the off diagonals.
[X,iter,del] = jaciter(A(P,:),B(P),1e-14,1000)
A\B
function [X,iter,del] = jaciter(A,b,tol,itmax)
% uses Jacobi iteration to solve the problem A*X==B.
%
% A: NxN array, presumed non-singular
% b: Nx1 vector
% tol: absolute tolerance for declaration of convergence
% itmax: maximum number of iterations allowed
%
% Note, the sizes of A and b are not checked for compatibility,
% nor is a check made for diagonal dominance.
n = size(A,1);
D = diag(A);
b = b(:);
LU = A - diag(D); % only the off-diagonal elements of A
Dinv = 1./D; % leave Dinv as a vector.
X0 = rand(n,1); % random start point for the iteration
del = inf;
iter = 0;
while (del>tol) && (iter <= itmax)
iter = iter + 1;
if iter >= itmax
warning('Maximum iterations exceeded')
end
% Note the use of .* below. This will cause an element-wise
% multiply to form X
X = Dinv.*(b - LU*X0);
del = norm(X - X0);
X0 = X;
end
if del > tol
warning('Convergence tolerance was not met at termination')
end
end
thank you for you're correction. Actually i get this matrix problem somewhre in a book called metode numerik by bambang triatmodjo page 85 (you can check it). 
I do find it strange when i see this matrix, because one of the conditions for using the Jacobi method is that the diagonal side of the matrix must be larger than the surrounding columns/rows. But I think that there is another solution to this problem. Once again, thank you for your solution. I do learn more from you, and i hope you can help me again in the future =)
Categories
Find more on Solver Outputs and Iterative Display in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
