The convergence of classical iterative solvers for linear systems is determined by the spectral radius of the iteration matrix, $\rho(\mathbf{G})$. For a general linear system, it is difficult to determine an optimal (or even good) SOR parameter due to the difficulty in determining the spectral radius of the iteration matrix. Below I have included many additional details, including an example of a real problem where the optimal SOR weight is known.
Spectral radius and convergence
The spectral radius is defined as the absolute value of the largest magnitude eigenvalue. A method will converge if $\rho<1$ and a smaller spectral radius means faster convergence. SOR works by altering the matrix splitting used to derive the iteration matrix based on the choice of a weighting parameter $\omega$, hopefully decreasing the spectral radius of the resulting iteration matrix.
Matrix splitting
For the discussion below, I will assume that the system to be solved is given by
$$\mathbf{A}x=b,$$
with an iteration of the form
$$x^{(k+1)}=v+\mathbf{G}x^{(k)},$$
where $v$ is a vector, and iteration number $k$ is denoted $x^{(k)}$.
SOR takes a weighted average of the old iteration and a Gauss-Seidel iteration. The Gauss-Seidel method relies on a matrix splitting of the form
$$\mathbf{A} = \mathbf{D} + \mathbf{L} + \mathbf{U}$$
where $\mathbf{D}$ is the diagonal of $\mathbf{A}$, $\mathbf{L}$ is a lower triangular matrix containing all elements of $\mathbf{A}$ strictly below the diagonal and $\mathbf{R}$ is an upper triangular matrix containing all elements of $\mathbf{A}$ strictly above the diagonal. The Gauss-Seidel iteration is then given by
$$x^{(k+1)}=(\mathbf{D}+\mathbf{L})^{-1}b+\mathbf{G}_{\rm G-S}x^{(k)}$$
and the iteration matrix is
$$\mathbf{G}_{\rm G-S}=-(\mathbf{D}+\mathbf{L})^{-1}\mathbf{U}.$$
SOR can then be written as
$$x^{(k+1)}=\omega(\mathbf{D}+\omega\mathbf{L})^{-1}b+\mathbf{G}_{\rm SOR}x^{(k)}$$
where
$$\mathbf{G}_{\rm SOR} = (\mathbf{D}+\omega\mathbf{L})^{-1}((1-\omega)\mathbf{D}-\omega\mathbf{U}).$$
Determining the convergence rate of the iterative scheme really boils down to determining the spectral radius of these iteration matrices. In general, this is a hard problem unless you know something specific about the structure of the matrix. There are very few examples that I'm aware of where the optimal weighting coefficient is computable. In practice, $\omega$ must be determined on the fly based on the observed (presumed) convergence of the running algorithm. This works in some cases, but fails in others.
Optimal SOR
One realistic example where the optimal weighting coefficient is known arises in the context of solving a Poisson equation:
$$\nabla^2 u = f~{\rm in}~\Omega\\
u = g~{\rm on}~\partial\Omega$$
Discretizing this system on a square domain in 2D using second order finite differences with uniform grid spacing results in a symmetric banded matrix with 4 on the diagonal, -1 immediately above and below the diagonal, and two more bands of -1 some distance from the diagonal. There are some differences due to boundary conditions, but that is the basic structure. Given this matrix, the provably optimal choice for SOR coefficient is given by
$$\omega = \dfrac{2}{1+\sin(\pi \Delta x/L)}$$
where $\Delta x$ is the grid spacing and $L$ is the domain size. Doing so for a simple case with a known solution gives the following error versus iteration number for these two methods:
As you can see, SOR reaches machine precision in about 100 iterations at which point Gauss-Seidel is about 25 orders of magnitude worse. If you want to play around with this example, I've included the MATLAB code I used below.
clear all
close all
%number of iterations:
niter = 150;
%number of grid points in each direction
N = 16;
% [x y] = ndgrid(linspace(0,1,N),linspace(0,1,N));
[x y] = ndgrid(linspace(-pi,pi,N),linspace(-pi,pi,N));
dx = x(2,1)-x(1,1);
L = x(N,1)-x(1,1);
%desired solution:
U = sin(x/2).*cos(y);
% Right hand side for the Poisson equation (computed from U to produce the
% desired known solution)
Ix = 2:N-1;
Iy = 2:N-1;
f = zeros(size(U));
f(Ix,Iy) = (-4*U(Ix,Iy)+U(Ix-1,Iy)+U(Ix+1,Iy)+U(Ix,Iy-1)+U(Ix,Iy+1));
figure(1)
clf
contourf(x,y,U,50,'linestyle','none')
title('True solution')
%initial guess (must match boundary conditions)
U0 = U;
U0(Ix,Iy) = rand(N-2);
%Gauss-Seidel iteration:
UGS = U0; EGS = zeros(1,niter);
for iter=1:niter
for iy=2:N-1
for ix=2:N-1
UGS(ix,iy) = -1/4*(f(ix,iy)-UGS(ix-1,iy)-UGS(ix+1,iy)-UGS(ix,iy-1)-UGS(ix,iy+1));
end
end
%error:
EGS(iter) = sum(sum((U-UGS).^2))/sum(sum(U.^2));
end
figure(2)
clf
contourf(x,y,UGS,50,'linestyle','none')
title(sprintf('Gauss-Seidel approximate solution, iteration %d', iter))
drawnow
%SOR iteration:
USOR = U0; ESOR = zeros(1,niter);
w = 2/(1+sin(pi*dx/L));
for iter=1:niter
for iy=2:N-1
for ix=2:N-1
USOR(ix,iy) = (1-w)*USOR(ix,iy)-w/4*(f(ix,iy)-USOR(ix-1,iy)-USOR(ix+1,iy)-USOR(ix,iy-1)-USOR(ix,iy+1));
end
end
%error:
ESOR(iter) = sum(sum((U-USOR).^2))/sum(sum(U.^2));
end
figure(4)
clf
contourf(x,y,USOR,50,'linestyle','none')
title(sprintf('Gauss-Seidel approximate solution, iteration %d', iter))
drawnow
figure(5)
clf
semilogy(EGS,'b')
hold on
semilogy(ESOR,'r')
title('L2 relative error')
xlabel('Iteration number')
legend('Gauss-Seidel','SOR','location','southwest')
