function [ errjac, errgs, errsor ] = liniter ( N, nit, iplot );
% LinIter
%
% This m-file is part of the material of the lecture
%   Numerical Mathematics
% at the Math Department of Halle University, Germany.
% Winter term 2004/05
%
% Author :      Prof. Dr. M. Arnold, arnold@mathematik.uni-halle.de
% Version of :  Nov 29, 2004
%
% Visualization of iterative methods for the solution of linear systems
%
% Parameters:
%   N      (input)  : number of grid points in FD approximation
%   nit    (input)  : number of iteration steps
%   iplot  (input)  : control flag "graphical output"
%                       =0 .. no graphical output
%                       =1 .. 3D animation
%   errjac (output) : error of Jacobi iteration (discrete l_2 norm)
%   errgs  (output) : error of Gauss-Seidel iteration (discrete l_2 norm)
%   errsor (output) : error of SOR iteration with optimal \omega (discrete l_2 norm)
%
% Examples:
%
%   [ errjac, errgs, errsor ] = liniter ( 20, 250, 1 );
%
%   see also plotiter.m
%
% Reference: Section 8.4 of
%   Stoer, J.; Bulirsch, R.: Numerische Mathematik 2, Dritte Auflage.
%   Springer-Verlag Berlin Heidelberg New York, 1990.

% -> finite difference approximation of Poisson problem in 2D
h = 1/N;
x = linspace ( 0, 1, N+1 );
[ XX, YY ] = meshgrid ( x );

% -> exact solution
e  = ones ( (N-1)^2, 1 );
ee = ones ( (N-1)^2, 1 );
ee(N-1:N-1:(N-1)^2) = 0;
L = spdiags ( [ e ee ] / h^2, [ -(N-1) -1 ], (N-1)^2, (N-1)^2 );
A = L + spdiags ( -4*e / h^2, 0, (N-1)^2, (N-1)^2 ) + L';

xsol = A \ ones( (N-1)^2, 1 );
Usol = zeros ( N+1, N+1 );
for i=1:N-1,
    Usol(1+(1:N-1),1+i) = xsol((i-1)*(N-1)+(1:N-1));
end;

umax = max ( max ( -Usol ) );

% -> initialize iteration
errjac = zeros ( nit+1, 1 );
errgs  = zeros ( nit+1, 1 );
errsor = zeros ( nit+1, 1 );

Ujac = zeros ( N+1, N+1 );
Ugs  = zeros ( N+1, N+1 );
Usor = zeros ( N+1, N+1 );

omega = 2 / ( 1 + sin ( pi / N ) );

% -> iteration
for it=0:nit,
    
    % -> evaluate error of actual iterates
    errjac(1+it) = sqrt ( sum ( sum ( ( Ujac-Usol ).^2 ) ) / (N-1)^2 ); 
    errgs(1+it)  = sqrt ( sum ( sum ( ( Ugs -Usol ).^2 ) ) / (N-1)^2 ); 
    errsor(1+it) = sqrt ( sum ( sum ( ( Usor-Usol ).^2 ) ) / (N-1)^2 ); 
    
    % -> visualization (if applicable)
    if iplot==1,
        for isub=1:4,
            subplot ( 2, 2, isub );  
            if isub==1,
                surf ( XX, YY, -Ujac );  title ( 'Jacobi iteration' );
            elseif isub==2,
                surf ( XX, YY, -Ugs  );  title ( 'Gauss-Seidel iteration' );
            elseif isub==3,
                surf ( XX, YY, -Usor );  title ( 'SOR method' ); 
            elseif isub==4,
                surf ( XX, YY, -Usol );  title ( 'Exact solution' );
            end;    
            xlabel ( 'x' );  ylabel ( 'y' );  zlabel ( '- u(x,y)' );
            set ( gca, 'ZLim', [ 0, 1.0*umax ] );
            set ( gca, 'CLim', [ 0,     umax ] );
        end;
        set ( gcf, 'Name', sprintf ( 'Iteration step   k = %i', it ) );
        drawnow,
    end;
    
    % -> iteration (max. nit iteration steps)
    if it<nit,
        
        % -> Jacobi
        Unew = zeros ( N+1, N+1 );
        for i=1:N-1,  for j=1:N-1,
            Unew(1+i,1+j) = ( h^2                                       ...
                            - Ujac(1+(i-1),1+j) - Ujac(1+(i+1),1+j)     ...
                            - Ujac(1+i,1+(j-1)) - Ujac(1+i,1+(j+1)) ) / (-4);
        end;  end;
        Ujac = Unew;
        
        % -> Gauss-Seidel
        for i=1:N-1,  for j=1:N-1,
            Ugs(1+i,1+j) = ( h^2                                     ...
                           - Ugs(1+(i-1),1+j) - Ugs(1+(i+1),1+j)     ...
                           - Ugs(1+i,1+(j-1)) - Ugs(1+i,1+(j+1)) ) / (-4);
        end;  end;
        
        % -> SOR
        for i=1:N-1,  for j=1:N-1,
            Usor(1+i,1+j) = ( 1 - omega ) * Usor(1+i,1+j)              ...
                           + omega * ( h^2                             ...
                           - Usor(1+(i-1),1+j) - Usor(1+(i+1),1+j)     ...
                           - Usor(1+i,1+(j-1)) - Usor(1+i,1+(j+1)) ) / (-4);
        end;  end;
        
    end;
    
end;    
