1

I need to code the Gauss Seidel and Successive over relaxation iterative methods in Matlab. I have created the below code for each of them, however my final solution vector does not return the correct answers and i'm really struggling to figure out why. Could anyone please help me? In both cases, x is the final solution vector and i returns the number of iterations.

Thanks in advance

Gauss Seidel Method:

function [x,i] = gaussSeidel(A,b,x0,tol)
x2 = x0;
count = 0;
D = diag(diag(A));
U = triu(A-D);
disp(U);
L = tril(A-D);
disp(L);
C = diag(diag(A));
disp(C);
Inv = inv(C+D);
error = inf;
      while error>tol
          x1 = x2;
          x2 = Inv*(b-(U*x1));
          error = max(abs(x2-x1)/abs(x1));
          count = count + 1;
      end
x = x2;
i = count;
end

SOR Method:

function [x,i] = sor(A,b,x0,tol,omega)
[m,n] = size(A);
D = diag(diag(A));
U = triu(A-D);
L = tril(A-D);
count = 1;
xtable = x0;
w = omega;
if size(b) ~= size(x0)
    error('The given approximation vector does not match the x vector size');
elseif m~=n
    error('The given coefficient matrix is not a square');
else
    xnew = (inv(D+w*L))*(((1-w)*D-w*U)*x0 +w*b);
    RelError = (abs(xnew-x0))/(abs(xnew));
    RelErrorCol = max(max(RelError));
    while RelErrorCol>tol
        xnew = (inv(D+w*L))*(((1-w)*D-w*U)*x0 +w*b);
        RelError = (abs(xnew-x0))/(abs(xnew));
        RelErrorCol = max(max(RelError));
        x0 = xnew;
        count = count+1;
        xtable = [xtable, xnew];
    end
    disp(xtable);
    x = xnew;
    i = count;
end
6
  • In your Gauss--Seidel function, there is a mistake: C and D are both equal to a diagonal matrix whose diagonal is that of A. That results in Inv being the inverse of 2*diag(diag(A)). According to the (standard) Gauss--Seidel algorithm, your Inv should be the inverse of A-U, where U is the matrix you compute.
    – Drake
    Apr 4, 2014 at 7:17
  • You may also want to consult this MatLab code from the Wikipedia entry on Gauss--Seidel
    – Drake
    Apr 4, 2014 at 7:40
  • Btw, in case you want to understand and fix your own code, the way you calculate the error is not correct: you end up taking the max of a matrix, which results in a vector used to check the precondition of the while loop.
    – Drake
    Apr 4, 2014 at 7:47
  • @Drake: Marko wants to solve Gauss-Seidel using matrix algebra instead of the point-based method. Good reference though.
    – rayryeng
    Apr 4, 2014 at 16:48
  • @rayryeng Where is this mentioned in the question? The only constraint I see is "code [...] in matlab".
    – Drake
    Apr 4, 2014 at 16:56

1 Answer 1

3

Gauss-Seidel: Your line that describes C is wrong. Actually it shouldn't be there. Also for the Inv line, it should be inv(D+L), not inv(C+D).

As for the SOR method, in hindsight it seems right. To double check, compare with this method:

http://www.netlib.org/templates/matlab/sor.m. This method relies on http://www.netlib.org/templates/matlab/split.m


Edit: April 4, 2014 - Also check: https://www.dropbox.com/s/p9wlzi9x9evqj5k/MTH719W2013_Assn4_Part1.pdf?dl=1 . I taught a course on Applied Linear Algebra and have MATLAB code that implements Gauss-Seidel and SOR. Check slides 12-20 for the theory and how to implement Gauss-Seidel and slides 35-37 for the SOR method.

Let me know how it goes.

4

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.