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
Marko Vidalis
  • 331
  • 1
  • 7
  • 18
  • 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 04 '14 at 07:17
  • You may also want to consult this [MatLab code from the Wikipedia entry on Gauss--Seidel](http://en.wikipedia.org/wiki/Gauss–Seidel_method#Program_to_solve_arbitrary_no._of_equations_using_Matlab) – Drake Apr 04 '14 at 07: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 04 '14 at 07:47
  • @Drake: Marko wants to solve Gauss-Seidel using matrix algebra instead of the point-based method. Good reference though. – rayryeng Apr 04 '14 at 16:48
  • @rayryeng Where is this mentioned in the question? The only constraint I see is "code [...] in matlab". – Drake Apr 04 '14 at 16:56
  • Sorry, I should have been more clear. In the code itself it is iteratively solving for the solution using matrix algebra. Using the upper and lower triangular matrices of the off-diagonals is one step in the solution. Suggest that the question be edited so that the desired method to be used for Gauss-Seidel and SOR is by using matrix algebra instead of point-based methods.... unless the person who asked the question can solve using any method. – rayryeng Apr 04 '14 at 17:03

1 Answers1

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.

rayryeng
  • 102,964
  • 22
  • 184
  • 193
  • I made that change to my gauss Seidel code, but still with no luck. Can you see any other errors in my code? Thank you – Marko Vidalis Apr 04 '14 at 09:29
  • I tried the code above with http://pastebin.com/HG5DLXuh . I got errors, it never estimated the correct vector. Why? – Royi Jan 11 '15 at 16:54
  • @Drazick - http://stackoverflow.com/questions/24730993/jacobi-iteration-doesnt-end/24733027#24733027 - Read the commentary on why Gauss-Seidel and Jacobi may not work. Make sure your system is diagonally dominant or you can't use these iterative methods. – rayryeng Jan 11 '15 at 17:26
  • @rayryeng see my comment at: http://stackoverflow.com/questions/24730993/jacobi-iteration-doesnt-end/24733027#comment44182909_24733027 – Royi Jan 11 '15 at 19:18