Gauss seidel method not working for sparse matrices

I have been tasked with using the gauss seidel method for a large 14 x 14 sparse matrix. When I run my code I obtain a vector of NaN values. The user defined function I am using is as follows:

function [x] = gaussseidel(A, b, x0, maxiter, tolerance)

%gaussseidel solves sets of linear equations using the gauss seidel method

% Detailed explanation goes here

x = x0;

for iter = 1:maxiter

x_old = x;

% Iterate through each equation

for i = 1:length(b)

% Calculate the sum of other terms except the diagonal term

sigma = A(i, 1:i-1) * x(1:i-1) + A(i, i+1:end) * x_old(i+1:end);

% Update the current solution

x(i) = (b(i) - sigma) / A(i, i);

end

% Check for convergence

if abs(x - x_old) < tolerance

fprintf('Converged in %d iterations\n', iter);

break;

end

end

The following code was written to use the above function:

A = [-1,0,0,0,0,0,0,0,1,0,0,0,-1,-1;

0,-1,0,0,0,0,0,0,0,1,0,0,-1,0;

0,0,-1,0,0,0,0,0,0,0,1,0,1,-1;

0,0,0,-1,0,0,0,0,0,0,0,1,0,1;

1,0,0,0,-1,0,0,0,-1,0,0,0,0,0;

0,1,0,0,0,-1,0,0,0,-1,0,0,0,0;

0,0,1,0,0,0,-1,0,0,0,-1,0,0,0;

0,0,0,1,0,0,0,-1,0,0,0,-1,0,0;

0,0,0,0,0,0,0,0,-0.9,0,0,0,1,1;

0,0,0,0,0,1,0,0,0,-1,0,0,0,0;

0,0,0.85,0,0,0,0,0,0,0,-1,0,0,0;

0,0,0,-0.65,0,0,0,0,0,0,0,1,0,0;

0.1,0,0,0,-1,0,0,0,0,0,0,0,0,0;

0,0,0,0,0,0,0,0,0,0,0,0,-0.3,0.7];

% Right-hand side vector b

b = [-10; -20; 0; 0; 0; 0; 0; 0; 9; 0; 0; 0; 0; 0];

x0 = zeros(size(b));

maxiter = 10000;

tolerance = 1e-12;

x = gaussseidel(A,b,x0,maxiter,tolerance);

disp('Solution:')

disp(x)

Using the function and mmatrices A and b above the following answers should be obatined:

1.0989

26.1538

26.3736

8.4772

0.1099

13.0769

3.9560

2.9670

0.9890

13.0769

22.4176

5.5102

6.9231

2.9670

Would anyone have any suggestions on how my user defined function could be altereed ti generate the correct answers for the above matrix using the gauss seidel method.

Thank you!!!!

[Torsten]

A(13,13) is zero - thus you divide by zero. You should at least arrange (by row and/or column exchange) that the diagonal of A has no zero elements.

[John D'Errico]

A = [-1,0,0,0,0,0,0,0,1,0,0,0,-1,-1;
    0,-1,0,0,0,0,0,0,0,1,0,0,-1,0;
    0,0,-1,0,0,0,0,0,0,0,1,0,1,-1;
    0,0,0,-1,0,0,0,0,0,0,0,1,0,1;
    1,0,0,0,-1,0,0,0,-1,0,0,0,0,0;
    0,1,0,0,0,-1,0,0,0,-1,0,0,0,0;
    0,0,1,0,0,0,-1,0,0,0,-1,0,0,0;
    0,0,0,1,0,0,0,-1,0,0,0,-1,0,0;
    0,0,0,0,0,0,0,0,-0.9,0,0,0,1,1;
    0,0,0,0,0,1,0,0,0,-1,0,0,0,0;
    0,0,0.85,0,0,0,0,0,0,0,-1,0,0,0;
    0,0,0,-0.65,0,0,0,0,0,0,0,1,0,0;
    0.1,0,0,0,-1,0,0,0,0,0,0,0,0,0;
    0,0,0,0,0,0,0,0,0,0,0,0,-0.3,0.7];
spy(A)

A is not singular. However, as Torsten points out, it has a zero pivot in row 13. So you could exchange row 13 with row 1, for example. (Don't forget to exchange the same rows of b too, if you do.)

HOWEVER, even if you do exactly that, I will still expect Gauss-Seiel to diverge on the row-permuted system. (I actually tested that fact out, and it does diverge.) In fact, if you do exchange rows 1 and 13, now the problem will be strongly not diagonally dominant, since element 1 of row 13 is tiny.

Answer 1

[John D'Errico]

Um, don't use Gauss-Seidel?

It is not that this matrix is sparse. (In fact, it was not even created as a sparse matrix. Sparsity has nothing to do with anything here, and not even big or sparse enough to be worth making into a sparse matrix.) It is that the matrix has properties that will cause it to be divergent for Gauss-Seidel.

It would be my guess that you were expected to see divergence for Gauss-Seidel, since this is almost certainly a homework assignment. Teachers use assignments like that as a teaching opportunity.

Note that your matrix is not diagonally dominant. Not even that close. I would expect to see divergence. Really.

Anyway, if this is a real problem that you want to solve for some other purpose than explicitly homework, then why not use backslash, the tool you SHOULD be using in the first place?

A = [-1,0,0,0,0,0,0,0,1,0,0,0,-1,-1;
    0,-1,0,0,0,0,0,0,0,1,0,0,-1,0;
    0,0,-1,0,0,0,0,0,0,0,1,0,1,-1;
    0,0,0,-1,0,0,0,0,0,0,0,1,0,1;
    1,0,0,0,-1,0,0,0,-1,0,0,0,0,0;
    0,1,0,0,0,-1,0,0,0,-1,0,0,0,0;
    0,0,1,0,0,0,-1,0,0,0,-1,0,0,0;
    0,0,0,1,0,0,0,-1,0,0,0,-1,0,0;
    0,0,0,0,0,0,0,0,-0.9,0,0,0,1,1;
    0,0,0,0,0,1,0,0,0,-1,0,0,0,0;
    0,0,0.85,0,0,0,0,0,0,0,-1,0,0,0;
    0,0,0,-0.65,0,0,0,0,0,0,0,1,0,0;
    0.1,0,0,0,-1,0,0,0,0,0,0,0,0,0;
    0,0,0,0,0,0,0,0,0,0,0,0,-0.3,0.7];
% Right-hand side vector b
b = [-10; -20; 0; 0; 0; 0; 0; 0; 9; 0; 0; 0; 0; 0];

x = A\b

x = 14×1

1.0989 26.1538 26.3736 8.4772 0.1099 13.0769 3.9560 2.9670 0.9890 13.0769

https://en.wikipedia.org/wiki/Gauss–Seidel_method

Is your matrix symmetric and positive definite?

eig(A)

ans = 
  -1.0000 + 0.0000i
  -1.0000 + 0.0000i
  -0.0780 + 0.0000i
  -1.9220 + 0.0000i
  -1.5200 + 1.1088i
  -1.5200 - 1.1088i
  -1.3412 + 1.1615i
  -1.3412 - 1.1615i
   0.6631 + 0.2982i
   0.6631 - 0.2982i
   0.5916 + 0.0000i
  -0.3177 + 0.0000i
  -0.4863 + 0.0000i
  -0.5916 + 0.0000i

As I said, certainly not any of the things that would create a convergent solve under Gauss-Seidel. I won't even spend the time to write my own Gauss-Seidel to verify divergence. (Ok, actually I did write a GS iteration code, since this question does get asked often enough. Even with row permutations to generate a problem with non-zero pivots on the diagonal, it still diverges.)

[John]

Hi John, thanks for your response. For this assignment we were explicitly asked to use the gauss seidel method

[John D'Errico]

Fine. But that does NOT mean you will expect to see it converge! As I said, it is often the case that you will be given a problem that is expected to diverge. Think of it as a learning experience, to learn what happens when you use the wrong method.

Another common case where teachers seem to do this (not me, I'd never do that to some poor, unsuspecting student) is to give them an ODE to solve, where they are EXPLICITLY told to use Euler's method. Since it is trivial to write a problem where Euler will fail miserably, this is an easy thing to do.

Sorry. You cannot force Gauss-Seidel to converge on such a problem. Ok, one thing I have never tried is to pray to my computer. I suppose it cannot hurt. So maybe try it. Be careful though, as you don't want to give your computer the impression that it is god. These days, with AI, that is probably a bad thing to do.