Unable to find solution to matrix using Gauss Seidal code. How should I proceed to get the solution?

4 views (last 30 days)
%Below is my program for Gauss Seidal method
A=[-0.707 -1 0 0 0 0 0 -1 0 0;
-0.707 0 0 0 0 0 0 0 -1 0;
0.707 0 0.707 0 0 0 0 0 0 0;
0.707 0 -0.707 -1 0 0 0 0 0 0;
0 0 0 1 0.5 -0.866 0 0 0 0;
0 0 0 0 0.866 0.5 0 0 0 0;
0 0 0 0 0 -0.866 -1 0 0 0;
0 0 0 0 0 -0.5 0 0 0 -1;
0 1 0.707 0 -0.5 0 -1 0 0 0;
0 0 -0.707 0 -0.866 0 0 0 0 0];
b=[0;0;-400;0;0;-200;0;0;0;0];
%Gauss Seidel
x0=[0;0;0;0;0;0;0;0;0;0];
max=10000;
tol=1e-4;
[n,m]=size(A);
x=x0;
C=-A;
for i=1:n
C(i,i)=0;
end
for i=1:n
C(i,:)=C(i,:)/A(i,i);
end
for i=1:n
d(i,1)=b(i)/A(i,i);
end
i=1;
while i<=max
xold=x;
for j=1:n
x(j) = C(j,:)*x + d(j);
end
if norm(x-xold) <= tol
disp('Gauss-Siedel method converge');
disp([x'])
return
end
i = i + 1;
end
disp('Method did not Converge')
Solution_by_Gauss_Seidal=x

Accepted Answer

John D'Errico
John D'Errico on 29 Sep 2021
Edited: John D'Errico on 29 Sep 2021
By the way, the name is spelled Seidel, not Seidal.
A=[-0.707 -1 0 0 0 0 0 -1 0 0;
-0.707 0 0 0 0 0 0 0 -1 0;
0.707 0 0.707 0 0 0 0 0 0 0;
0.707 0 -0.707 -1 0 0 0 0 0 0;
0 0 0 1 0.5 -0.866 0 0 0 0;
0 0 0 0 0.866 0.5 0 0 0 0;
0 0 0 0 0 -0.866 -1 0 0 0;
0 0 0 0 0 -0.5 0 0 0 -1;
0 1 0.707 0 -0.5 0 -1 0 0 0;
0 0 -0.707 0 -0.866 0 0 0 0 0];
First, is A singular?
cond(A)
ans = 7.4663
So it is quite well conditioned, not even close to singular. But does that mean Gauss-Seidel will always converge? Now look at the same page about Gauss-Seidel. Do some reading.
In there, we will see the comment:
The convergence properties of the Gauss–Seidel method are dependent on the matrix A. Namely, the procedure is known to converge if either:
The Gauss–Seidel method sometimes converges even if these conditions are not satisfied.
However, those zeros on the diagonal are complete killers.
Is A symmetric positive definite? NO. It is not symmetric at all, and it is cerrainly not positive definite. Is A even diagonally dominant? NO. Again, if necessary, do some reading. What does diagonally dominat mean?
The zeros on the diagonal kill you.
diag(A)
ans = 10×1
-0.7070 0 0.7070 -1.0000 0.5000 0.5000 -1.0000 0 0 0
Not all matrix problems can be solved using Gauss-Seidel. Maybe Gauss-Seidal will work though. :)
Does that mean the problem cannot be solved? Of course not. A is non-singular. The solution is trivial.
b=[0;0;-400;0;0;-200;0;0;0;0];
X = A\b
X = 10×1
-548.1782 387.5620 -17.5927 -375.1240 14.3626 -424.8760 367.9427 0.0000 387.5620 212.4380
Can you solve the problem using gauss-Seidel? Well, not directly, no. Can you solve it IF you permute the columns of A? For example...
p = [8 9 3 1 4 6 7 10 2 5];
spy(A(:,p))
diag(A)
ans = 10×1
-0.7070 0 0.7070 -1.0000 0.5000 0.5000 -1.0000 0 0 0
This column reordering implies an implicit reordering of the unknown elements of x. Will it be convergent for Gauss-Seidel? Possibly. Then you would need to re-order the unknowns.
Or, you could try to find a row-re-ordering of A, such that A is again non-zero on the diagonal. Then you would reorder b in the same way. That will NOT require permuting the unknowns at the end.
r = [3 9 10 4 5 6 7 1 2 8];
spy(A(r,:))
diag(A(r,:))
ans = 10×1
0.7070 1.0000 -0.7070 -1.0000 0.5000 0.5000 -1.0000 -1.0000 -1.0000 -1.0000
We can see that
A(r,:)\b(r)
ans = 10×1
-548.1782 387.5620 -17.5927 -375.1240 14.3626 -424.8760 367.9427 -0.0000 387.5620 212.4380
yields the same solution as A\b. Will this be convergent for your Gauss-Seidel code? Perhaps. At least it may have a chance of converging. Will it? Possibly.
What is the spectral radius of A(r,:)? It needs to be less than 1 for Gauss-seidel to have assured convergence. Here, it seems to be a bit larger than 1.
format long g
max(abs(eig(A(r,:))))
ans =
1.15399258885417
Not too far off, but you may still have problems.
In fact, there are row permutations of A that will reduce the spectral radius of A to just exactly 1.
Is there a good reason why you think you NEED to use Gauss-Seidel?
  1 Comment
John D'Errico
John D'Errico on 1 Oct 2021
@Nilotpal Kalita: Please don't use an answer to make a comment. You said:
Thank you for the detailed explanation.
I was trying to compare this method with a Gauss elimination code."
My response is not all linear systems are solvable using Gauss-Seidel. Diagonal dominace is useful in this respect. The typical limit is that the spectral radius must be no greater than 1.
A=[-0.707 -1 0 0 0 0 0 -1 0 0;
-0.707 0 0 0 0 0 0 0 -1 0;
0.707 0 0.707 0 0 0 0 0 0 0;
0.707 0 -0.707 -1 0 0 0 0 0 0;
0 0 0 1 0.5 -0.866 0 0 0 0;
0 0 0 0 0.866 0.5 0 0 0 0;
0 0 0 0 0 -0.866 -1 0 0 0;
0 0 0 0 0 -0.5 0 0 0 -1;
0 1 0.707 0 -0.5 0 -1 0 0 0;
0 0 -0.707 0 -0.866 0 0 0 0 0];
max(abs(eig(A)))
ans = 1.2997
Which clearly exceeds 1. And that tells me that no matter what, a code like Gauss-Seidel will not be successful with this matrix as it is. I did show, that with some effort, you can permute the rows or the columns of A to enable success. But seriously, why bother? Gauss-Seidel is never a good solution method to solve a linear system of equations. Yes, they teach it in school, essentially as a way to talk about the solution of linear systems of equations. The general idea can even be useful for nonlinear problems. But you should not be seriously thinking about using it. For that matter, you should not even be writing your own gaussian elimination code. If you need to solve a linear system, there are multiple ways provided in MATLAB, written by professionals, who know what they are doing.
Start with backslash, or LINSOLVE. If you need iterative schemes, then look at tools like LSQR, etc. If you need to use factorization schemes, use tools provided like LU or QR.

Sign in to comment.

More Answers (2)

Alan Stevens
Alan Stevens on 29 Sep 2021
You are dividing by A(i,i) some of which are zero. These will introduce NaNs.

Nilotpal Kalita
Nilotpal Kalita on 30 Sep 2021
Thank you for the detailed explanation.
I was trying to compare this method with a Gauss elimination code.

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!