Implementing Gauss Seidel iterative solver

Here I would use notation Ax=b for system of linear equations, where A is a square matrix, x is vector of unknowns and b is vector of right hand side values.

Gaussian elimination as direct method for solving large system of linear equations would take 2*(n^3/3) + O(n^2) time, so it has roughly O(n^3) complexity.
Cholesky decomposition is only two times faster than Gaussian elimination so it is also O(n^3) while A must be positive definite for this method.

Thus for large systems only iterative methods are used.

As such there is Gauss-Seidel iterative process which converges to real solution in cases when A is diagonally dominant. For example, such matrix occurs as result of finite difference approximation for solving PDEs (see Poisson equation). The iterative process can be stopped after some finite time when it gives result close to the exact solution.
For checking how close is the current solution to the exact values one can use residual (res=b-Ax) norms.

Gauss Seidel method can be defined in two ways although there are equivalent.

In books as well as on the web, there are explicit formula in matrix terms and also explicit formula expressing components of x.

Also see graphical representation of Gauss-Seidel iterations.

Here is Matlab code for Gauss Seidel in matrix terms and also Matlab file for Gauss Seidel as a cycle for x values. There is also example code which solves two dimensional stationary heat equation.

It turns out is that matrix version is two times faster. This is because Matlab is optimized to perform matrix operations efficiently.

But here is important thing to mention.

When there is a sparse matrix, one should consider implementing tailored code which avoids storing coefficients matrix in memory. For example, PDE problems create discretization matrix which has 3-5 or other limited number of different values (depending on discretization stencil).
These values can be hardcoded in the implementation and this gives you hundred times better performance then using same algorithm with general input requirements.

One Response to “Implementing Gauss Seidel iterative solver”

  1. Anie Ekpes says:

    Please when i tried ruuning the “Gauss seidel in matrix terms” in MATLAB, I got this error message:

    ??? Input argument “A” is undefined.

    Error in ==> gseidel1 at 13
    n = size(A,1);

    Please how do I get rid of this Error/Debug it?
    Thanks, please email me your response, I am Using MATLAB 7b

RSS feed for comments on this post. And trackBack URL.

Leave a Reply