i was calculate using gauss jacobi method with 4x4 matrice, why does NaN appear? what should i do?

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]);

Answer 1

[Torsten]

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)

abcd
         0         0         0         0
    0.5000    1.0000    1.5000   -3.2000
    1.3500   -0.7333    2.5500   -2.7000
    0.5750   -1.2000    1.9917   -2.8367
    0.9225   -0.8011    1.9092   -3.0417
    1.0663   -0.9578    2.0601   -2.9784
    0.9591   -1.0349    2.0052   -2.9795
    0.9872   -0.9813    1.9862   -3.0060
    1.0099   -0.9931    2.0062   -2.9990
    0.9964   -1.0050    2.0015   -2.9974
    0.9980   -0.9984    1.9981   -3.0007
    1.0013   -0.9989    2.0006   -3.0001
    0.9997   -1.0006    2.0003   -2.9997
    0.9997   -0.9999    1.9998   -3.0001
    1.0002   -0.9998    2.0000   -3.0000
    1.0000   -1.0001    2.0000   -3.0000
    1.0000   -1.0000    2.0000   -3.0000
    1.0000   -1.0000    2.0000   -3.0000
    1.0000   -1.0000    2.0000   -3.0000
    1.0000   -1.0000    2.0000   -3.0000
    1.0000   -1.0000    2.0000   -3.0000

eror
   2.3864e-05

Output argument "H" (and possibly others) not assigned a value in the execution with "solution>jacobi" function.

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

[Beatrice]

thank you for your solution. I can learn more from you, and i hope you can help me again in the future =)

Answer 2

[John D'Errico]

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]

A = 4×4

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.

https://en.wikipedia.org/wiki/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.

[Torsten]

@John D'Errico

After row permutation, the Jacobi method can be applied.

[John D'Errico]

@Torsten

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)'

A = 3×3

1 2 3 4 5 6 7 8 9

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)

ans = 2

A(3,3) = 9 - 100*eps(9);
rank(A)

ans = 3

Even a broader change will suffice for failure of the Jacobi method.

A(3,3) = 8;
rank(A)

ans = 3

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)

A2 = 3×3

11 2 3 4 15 6 7 8 18

b = [2;3;5];
A2\b % known solution

ans = 3×1

0.1123 0.0929 0.1928

[X,iter,del] = jaciter(A2,b,1e-14,1000)

X = 3×1

0.1123 0.0929 0.1928

iter = 78

del = 7.7778e-15

So it converged to the known solution, after 76 iterations. Yet on problem A, it diverges

A\b % known solution

ans = 3×1

-2.3333 3.6667 -1.0000

[X,iter,del] = jaciter(A,b,1e-14,1000)

X = 3×1

NaN NaN -Inf

iter = 723

del = NaN

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],:)

A3 = 3×3

4 15 6 7 8 18 11 2 3

b3 = b([2 3 1])

b3 = 3×1

3 5 2

This does not change the system.

A3\b3

ans = 3×1

0.1123 0.0929 0.1928

Yet now jaciter should fail.

[X,iter,del] = jaciter(A3,b3,1e-14,1000)

X = 3×1

NaN NaN NaN

iter = 493

del = NaN

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]

After row permutation, the Jacobi method can be applied.

Well, I meant in the actual case, of course.

[John D'Errico]

@Torsten

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]

A = 4×4

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,:)

ans = 4×4

2 0 1 1 1 3 1 -1 0 -1 4 1 0 1 1 -5

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)

X = 4×1

1.0000 -1.0000 2.0000 -3.0000

iter = 50

del = 7.9813e-15

A\B

ans = 4×1

1.0000 -1.0000 2.0000 -3.0000

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

[Beatrice]

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 =)