Jacobi iterative method problem

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!

[Walter Roberson]

Please post sample a and b and optional arguments for us to test with.

[Agostino Dorano]

I paste here same example of use :

>> a = randi(5,5)

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 = randi(100,5,1)

b =

    76
    75
    40
    66
    18

>> a = sparse(a); % the function require a sparse matrix
>> jacobiSolver(a,b)
Warning: Limit reached, probably not convergence... calculation stopped and resrel = 0 
> In jacobiSolver (line 33) 

ans =

  1.0e+295 *

    1.7645
    5.8905
    3.5209
    4.7615
    5.0186
    
%% now i paste the result obteined with "\"
>> a\b

ans =

  101.4595
 -218.4595
 -286.4054
  515.7297
 -110.5405

Optional arguments of the function are NMAX and TOL both utilized in this instruction:

 while norm(x-x0,1)>=TOLX && nit<NMAX

but I don't use them in the example. Defaults value are 500 iteractions for NMAX and TOL = 10 ^(-6)

[Walter Roberson]

checkAB() ?

[Agostino Dorano]

checkAB(a,b) is a local function that check if a and b are correctly inserted by users. I paste here my code for check of All function arguments (all local functions):

function checkAB(a,b)
if(~ismatrix(a) ||  ~issparse(a) || length(size(a))~=2 || size(a,1)~=size(a,2) || size(a,1) < 1  )
    error("A must be a sparse square matrix with n greather than 1");
elseif (~isvector(b) ||~iscolumn(b) ||isa(a,'cell') || isa(a,'struct') || isa(a,'function_handle')|| isa(a,'table') || length(b)~=size(a,1)) 
    error("b must be a vector with dimension compatible with A");
elseif ~isnumeric(a) || ~isreal(a) || ~all(all(isfinite(a)))
    error("Elements of A must be real numbers different from + - inf or NaN");
elseif ~isnumeric(b) || ~isreal(b) || ~all(isfinite(b))
    error("Elements of b must be real numbers different from + - inf or NaN");
    %%% aggiungere condizioni tipiche di fallimento dei metodi iterativi
elseif ~(any(diag(a))) % da controllare => probabilmente mettere qualche norma
    error("null element on diag");
end
end

% CHECK FUNCTION OF NMAX PARAMTER
function [TOL, NMAX] = CheckNMAX(tol, nmax)
if (~isreal(nmax)||~isnumeric(nmax)||~(length(nmax)==1)||~isfinite(nmax)||nmax<=0 )
    error('NMAX must be a real number greather than 0.');
elseif ~(floor(nmax)==nmax)
    error("NMAX must be an integer value.");
end
[TOL,NMAX] = CheckTOL(tol,nmax);
end

% CHECK FUNCTION OF TOL PARAMETER
function [TOL, NMAX] = CheckTOL(tol,nmax)
if (~isreal(tol)||~isnumeric(tol)||~(length(tol)==1))
    error('TOL must be a real number.');
elseif ~isfinite(tol)
    error("TOL can't be + - Inf or NaN.");
elseif tol < eps || tol >=1 
    warning("tolerance value is wrong, see doc jacobiSolver. Setted to default value.");
    tol = 10^-6;
end
TOL = tol; NMAX = nmax;  
end

[Walter Roberson]

You calculate a mostly negative B and a positive c, and you have a recursion like

f(n) = B*f(n-1) + c

and you stop when norm(f(n)-f(n-1),1) < tolerance

but f(n) - f(n-1) is going to be B*f(n-1)+c - f(n-1) = (B-1)*f(n-1) + c and the 1 norm is the sum of the absolute value of the components, so sum( abs((B(K)-1)*x0(K) + c(K)) ) < tolerance . An x0 that would be true for would be x0 = -c./(B-1) or close by.

Anyhow, you start with an x0 and head for smaller and smaller x0 looking for one that will have the B vs c balance. However, your code with its 1-norm cannot tell which side of the balance you are on, so when you overshoot the balance, you just keep going in the same direction, getting further and further from the balance.

[valerio auricchio]

Sorry Walter Roberson can you explain bettere this "An x0 that would be true for would be x0 = -c./B or close by.

Anyhow, you start with an x0 and head for smaller and smaller x0 looking for one that will have the B vs c balance. However, your code with its 1-norm cannot tell which side of the balance you are on, so when you overshoot the balance, you just keep going in the same direction, getting further and further from the balance." I don't think I caught well.

[Walter Roberson]

If x0 exactly equaled -c./(B-1) then each abs((B(K)-1)*x0(K)+c(K)) component would be 0, and the sum would be 0. But you do not need an exact 0, so the x0 values can potentially be close to -c./(B-1) and still have the 1-norm less than the tolerance.

(Note: I made a correction for the -1)

Answer 1

[John D'Errico]

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.