MATLAB: Jacobi iterative method problem

equationiterationiterativejacobilinearmathematicsMATLABsolvesystem

I've implemented the Jacobi method in matlab but when i try it , the function give me wrongs results.. I don't know what I'm wrong. I post here my code and results.
function [x,niter,resrel] = jacobiSolver(a, b, varargin)
narginchk(2, 5); % check if the function receives the right number of input parameters
nargoutchk(0,3); % check if the function receives the right number of output parameters
checkAB(a,b) % priority control
% OPTIONAL PARAMETER CONTROL
switch nargin
    case 4; [TOL, NMAX] = CheckNMAX(varargin{1}, varargin{2});
    case 3; [TOL, NMAX] = CheckTOL(varargin{1}, 500);
    case 2; TOL = 10^-6; NMAX = 500;
    otherwise
end
% codifichiamo l'algoritmo => vedere se dobbiamo sostituire con spdiags
c = b./diag(a);
B = (-1./diag(a)')' .*(a-diag(diag(a))); 
x0=zeros(length(b),1);
x=c;
nit = 0;
%TOLX = max(TOL*norm(x,1),realmin);
TOLX = TOL;
while norm(x-x0,1)>=TOLX && nit<NMAX
    x0 = x ;
    x = B*x0 + c;
    TOLX = max(realmin, TOL*norm(x0,1));
    nit = nit+1
end
% IF THE NMAX VALUE IS REACHED THEN THE USER WILL BE NOTIFIED WITH A WARNING
if(nit==NMAX) 
    warning("Limit reached, probably not convergence... calculation stopped and resrel = "+0);
end
if nargout > 1
    niter = nit;
    if nargout == 3
        resrel = 0;
    end  
end
end
results obteined with simple matrix a and b:
Warning: Limit reached, probably not convergence... calculation stopped and resrel = 0 
> In jacobiSolver (line 32) 
ans =
  1.0e+295 *
    0.3636
    1.2137
    0.7255
    0.9811
    1.0341
correct results:
    ans =
   20.2162
  -43.2162
  -56.8378
  102.1081
  -21.7838
also are there better ways to implement this algorithm?
thanks a lot for helping!

Best Answer

Note that the Jacobi method for solution of a linear system is only convergent for SOME matrices. It is NOT always convergent. You might want to do some reading, here:
https://en.wikipedia.org/wiki/Jacobi_method
The point is, if your system is not diagonally dominant, then the Jacobi iteration will not converge. It will diverge. You are trying this on a RANDOM matrix. But a divergent result tends to imply divergence has occurred. Digging down into your comments, I see this matrix. Sigh. I wonder if it is diagonally dominant? Any takers on that bet?
a = [5     1     1     1     4
     5     2     5     3     1
     1     3     5     5     5
     5     5     3     4     5
     4     5     5     5     4];
b =[76
    75
    40
    66
    18];
So, what does diagonally dominant mean? The requirement is that
abs(A(i,i)) >= sum(abs(A(i,j)))
where that sum is taken over all j~=i.
Now, look at the 5th row of a. a(5,5) is 4. The other elements of a on that row are [4 5 5 5]. They sum to 19. Is 19>4? (Yes.) Therefore your matrix is NOT diagonally dominant. Therefore, you would expect the Jacobi iteration on this problem to at least potentially be divergent. In fact, none of the rows of a satisfy the requirement.
Did it diverge? Yes. Case closed. When you see a problem with a method diverging, don't just give up. Ask if it diverged for a valid reason. Do some reading about the method in question.
Suppose we tried a different matrix instead? I'll just use the guts of your code, because I don't feel like finding and saving all of your associated functions.
a = a + eye(5)*30;
c = b./diag(a);
B = (-1./diag(a)')' .*(a-diag(diag(a))); 
x0=zeros(length(b),1);
x=c;
nit = 0;
NMAX = 1000;
TOL = 1e-12;
TOLX = 1e-12;
while norm(x-x0,1)>=TOLX && nit<NMAX
    x0 = x ;
    x = B*x0 + c;
    TOLX = max(realmin, TOL*norm(x0,1));
    nit = nit+1;
end
Did it work?
x
x =
       2.0931
       1.7777
      0.77996
       1.3459
      -0.2909
>> a\b
ans =
       2.0931
       1.7777
      0.77996
       1.3459
      -0.2909
>> a
a =
    35     1     1     1     4
     5    32     5     3     1
     1     3    35     5     5
     5     5     3    34     5
     4     5     5     5    34
Of course. It worked because the Jacobi method is convergent for that matrix. Note that if you are taking linear algebra, and they are teaching the Jacobi iteration method for linear systems, this should be one of the things they would have dicussed at some point. I suppose it MIGHT be the thing they discuss AFTER you try using it to solve a non-convergent problem.
Related Question