echo on
%%%%% laplace_gs_sor.m %%%%%

disp(sprintf('**** Exampel of Gauss-Seidel iteration for Laplace eq. with SOR ****\n'));

%Two dimensional 3x3 elliptic boundary problem (Laplace equation)
%with boundary values 0 (left, bottom right) and 1 up.
%Natural ordering of nodes: 1 left bottom, ..., 3 right bottom,
%4 left second row (from bottom) etc. 9 right up.
%Laplace matrix L is composed of sub-matrices T and -I where
format short;
b=[0 0 0 0 0 0 1 1 1]'; %boundary conditions (upper edge = 1)
T=diag([4 4 4],0)+ diag([-1 -1],1)+ diag([-1 -1],-1)
I=eye(3)
Z=zeros(3,3)
LA=[T -I Z; -I T -I; Z -I T;]
%All diagonal elements are the same. Matrix LA can be saved by diagonals,
%or not at all because their structure is known.
%The system dimension (number of inner points, elements of vector x) is 9.
%Use an array 5x5 for the representation of the system domain with 
%initial and boundary conditions. 
%Ordering of nodes: (2,2) left up, ..., (2,4) right up,
% (3,2) left second row (from bottom) etc. (3,3) right bottom.
%Using Gauss-Seidel iteration no copy for results is needed:
%Initial approximation + boundary conditions
%D=[ones(1,5);0, [4 4 4], 0; 0, [4 4 4], 0;0, [4 4 4], 0;zeros(1,5);]
D=[ones(1,5);0, [0 0 0], 0; 0, [0 0 0], 0;0, [0 0 0], 0;zeros(1,5);]
U=D;	% copy for the old solution
om=1.1;	%relaxation parameter has to be determined
%flops(0);
%first approximation:
for i=2:4
  for j=2:4
     D(i,j)=U(i,j) + om*((D(i-1,j)+D(i,j-1)+D(i+1,j)+D(i,j+1))/4 - U(i,j))
  end
end
U=D;
%other steps
echo off;
for k=1:6	%number of additional steps
   for i=2:4
     for j=2:4
     	D(i,j)=U(i,j) + om*((D(i-1,j)+D(i,j-1)+D(i+1,j)+D(i,j+1))/4 - U(i,j));
     end
  end
  U=D;
  k=k+1;
end
echo on;
D	%approximate solution after 7 steps (faster convergence than 
		%Jacobi and Gauss-Seidel)
%flops
x=LA\b	%exact solution