This is a small and simple program to solve the linear system Au=F when A is a diagonally dominant matrix. Although the program is quite simple, I didn't find it in this page, so I decided to post it.
This program may be useful to people programing solving Partial Differential Equations via Finite Differences. This method may be more efficient than Gauss elimination for large matrices, usually found for 2D and 3D problems.
[u, it] = sor(A,F) finds the solution of the linear system applying successive
under/over relaxation technique. In this case the parameter w, the initial guess u0, the stopping criterion dtol and the maximum number of iterations
itmax are to default values: w = 1.5, u0 = F.*0+1, dtol = 1e-3, itmax = 1000.
The parameter w can theoretically take values between 0 - 2. However, the number of iterations may be reduced if is taken near 1.5. Using w = 1 will be equivalent to use Gauss-Seidel method.
u0 is the initial guess used by SOR to find the solution. The better the initial guess the faster the solution is found. However, if the matrix A is strictly diagonally dominant, the method converges to the solution no matter the initial guess (see
Lecture 6 "Iterative Techniques").
As stated before the matrix A should be diagonally dominant. This program DOESN'T check this property; you will just get a NaN vector or something similar as an answer.
Finally the program will stop when abs(u-ui)/abs(ui) < dtol, where u and
ui are the solution and initial guess used for a particular iteration.
I'd first comment that much of the time, the direct sparse solver in Matlab is a FAR better choice than this SOR code anyway.
n = 1000;
A = sprandn(n,n,.001)/10 + speye(n,n);
f = randn(n,1);
A is a sparse matrix which is diagonally dominant.
How fast is SOR compared to backslash?
So in this test case, backslash was 673 times as fast as SOR.
How about accuracy?
x0 = A\f;
xs = sor(A,f);
std(f - A*x0)
std(f - A*xs)
A difference of 12 extra digits seems important to me.
Next, a look at the help. It is decent, although not really in a standard form. I'd have preferred a more explicit set of comments that described each input argument, its shape, and its meaning, plus the default values.
More importantly, if you do
you will find NO explicit listing of the arguments, as they appear in the function header. You find these variables referenced in the help, but no place do you see in which order they must appear.
I did find an example of use, and an H1 line. There were reasonable error checks made (although not for the requirement of diagonal dominance, not hard to check) and default parameters were supplied.
One other comment - the suggestion is made that a value of 1.5 for the relaxation parameter is somehow optimal, or nearly so. As I recall, the optimal relaxation parameter is dependent on the system. A value of 1.5 may often be a reasonable choice.
The author might find the following code snippet of interest however. As applied to the simple test problem I used above, where do we see the best solution performance?
w = .5:.01:1.99;
err = zeros(size(w));
iter = err;
for i = 1:length(w)
[x,iter(i)] = sor(A,f,w(i));
err(i) = std(A*x-f);
Note that w = 1.5 did NOT yield the best result, neither in terms of speed or in accuracy. For the chosen problem, I found those to happen at roughly w = 1.