%
% apply the Jacobi, Gauss-Seidel and SOR method to compute an approximate 
% solution of the finite difference discretizatin of the BVP
%
%     - diff  y''(x) + adv y'(x) + react y(x) = f(x)  x in [a,b],   
%
%                          y(a) = r_a,
%                          y(b) = r_b,
%
% on [a,b] = [0,1].  The implementation is matrix-free, i.e., only the
% stencil is used.
%
% Matthias Heinkenschloss
% Jan 24, 2017   (revised Jan 22, 2020)
%

set(0, 'defaultaxesfontsize',24,'defaultaxeslinewidth',3.0,...
       'defaultlinelinewidth',3.0,'defaultpatchlinewidth',3.0,...
       'defaulttextfontsize',24);

clear all

%% Stopping tolerance for the iterative methods
tol = 1.e-4;

%% set problem data
example = 1;

% set mesh points
n   = 9;   % try n = 19, 99, 499 when upwind = 0, and n = 19, 49, 99 when upwind = 1.
hx  = 1/(n+1);
x   = (0:hx:1)';

if( example == 1 )
    % Laplace equation
    diff  = 1;
    adv   = 0;
    react = 0;
    ga    = 0;
    gb    = 0;
    
    f     = rand(size(x));
    yex   = zeros(size(x));
end
if( example == 2 )
    % advection diffusion equation
    diff  = 1.e-3;
    adv   = 1;
    react = 0;
    ga    = 0;
    gb    = 0;
    
    f     = ones(size(x));
    yex   = x - (exp((x-1)/diff)-exp(-1/diff))./(1 - exp(-1/diff));
end




%% Generate the 3-point stencil
% Stencil for Central Difference schme
%Stencil  = (1/hx^2)*[ -(diff + 0.5*hx*adv), (2*diff + hx*hx*react), -(diff - 0.5*hx*adv)]';
% Stencil for upwind scheme   
Stencil  = (1/hx^2)*[ -(diff + hx*adv), (2*diff + hx*adv + hx*hx*react), -diff]';
    
% stencil using in Jacobi and SOR iteration
Stenciliter = Stencil; Stenciliter(2) = 0;



% Set boundary data
y      = zeros(n+2,1);
y(1)   = ga;
y(n+2) = gb;

% starting value to use for all iterative methods
y0     = y; 
    
    
%% damped Jacobi method
omega = 1;
res   = zeros(n+2);
y     = y0;
iter  = 0;
ynew  = y;
for i = 2:n+1
    res(i) = sum( Stencil.*y(i-1:i+1) )  -  f(i);
end
normres = norm(res,2);
reshist = normres;
while normres > tol
    iter = iter+1;
    % use left-right ordering
    for i = 2:n+1
        ynew(i) = (1-omega)*y(i) ...
                  + omega*( f(i) - sum( Stenciliter.*y(i-1:i+1) ) )/ Stencil(2);
    end
    y = ynew;
    for i = 2:n+1
        res(i) = sum( Stencil.*y(i-1:i+1) )  -  f(i);
    end
    normres = norm(res(:),2);
    reshist = [ reshist; normres ];
    
%     figure(3)
%     plot(x,y,'-');
%     xlabel('x')
%     title(['Approximation in Iteration ',int2str(iter)])
%     disp('pause ...'); pause
end

figure(1)
semilogy(reshist,'d-'); hold on;
xlabel('Iteration')
ylabel('Residual Norm')



%% SOR method
% order nodes left to right
omega = 1;
y     = y0;
iter  = 0;
res   = zeros(n+2);
for i = 2:n+1
    res(i) = sum( Stencil.*y(i-1:i+1) )  -  f(i);
end
normres = norm(res(:),2);
reshist = normres;
while normres > tol
    iter = iter+1;
    % use lexicographic ordering
    for i = 2:n+1
        y(i) = (1-omega)*y(i) ...
               + omega*( f(i) - sum( Stenciliter.*y(i-1:i+1) ) )/ Stencil(2);
    end
    for i = 2:n+1
        res(i) = sum( Stencil.*y(i-1:i+1) )  -  f(i);
    end
    normres = norm(res(:),2);
    reshist = [ reshist; normres ];
    
%     figure(3)
%     plot(x,y,'-');
%     xlabel('x')
%     title(['Approximation in Iteration ',int2str(iter)])
%     disp('pause ...'); pause  
end

figure(1)
semilogy(reshist,'o-'); hold on;
xlabel('Iteration')
ylabel('Residual Norm')

% repeat using red-black ordering
y     = y0;
iter  = 0;
res   = zeros(n+2);
for i = 2:n+1
    res(i) = sum(sum(  Stencil.*y(i-1:i+1) ))  -  f(i);
end
normres = norm(res(:),2);
reshist = normres;
while normres > tol
    iter = iter+1;
    % use red-black
    for i = 2:2:n+1
        y(i) = (1-omega)*y(i) ...
               + omega*( f(i) - sum( Stenciliter.*y(i-1:i+1) ) )/ Stencil(2);
    end
    for i = 3:2:n+1
        y(i) = (1-omega)*y(i) ...
               + omega*( f(i) - sum( Stenciliter.*y(i-1:i+1) ) )/ Stencil(2);
    end
    for i = 2:n+1
        res(i) = sum( Stencil.*y(i-1:i+1) )  -  f(i);
    end
    normres = norm(res(:),2);
    reshist = [ reshist; normres ];
    
%     figure(3)
%     plot(x,y,'-');
%     xlabel('x')
%     title(['Approximation in Iteration ',int2str(iter)])
%     disp('pause ...'); pause
end

figure(1)
semilogy(reshist,'*--'); hold off;
xlabel('Iteration')
ylabel('Residual Norm')

legend('Jacobi','GS - lex', 'GS - rb')

figure(2)
plot(x,y,'-');
xlabel('x')
title('Computed solution ')

