% script file BdryLayer
% demonstrate the boundary layers 
% equation
% Du + D^2u = f
% u(a) = uA, u(b) = uB;
% and also analyse the properties of 
% linear system
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clc; clear all; %close all;

a = 0; b = 1;
n = 61;

lambda = -0.001; % diffusion

uA = 0; uB = 0;

x = linspace(a,b,n)';
h = x(2)-x(1);

main_diag = ones(n,1);

D = spdiags([-main_diag/h main_diag/h],[0 1],n,n);
D(1,1) = 0; D(1,2) = 0; 
D(n,n-1) = 0; D(n,n) = 0;

A = spdiags([main_diag/h/h -2*main_diag/h/h main_diag/h/h],[-1 0 1],n,n);
A(1,1) = -1; A(1,2) = 0;
A(n,n-1) = 0; A(n,n) = -1;

K = D+lambda*A;

% right hand side at the grid points
f = -ones(size(x));
% at the boundary points overwrite values
f(1) = -uA; f(n) = -uB;

u = K\f;

figure(1); clf;
subplot(1,2,1);
plot(x,u,'.-');
%axis equal

E = eig(full(K));
subplot(1,2,2);
plot(real(E),imag(E),'.');
%axis equal

%return

% Jacobi iteration
maxit = 1000;
toler = 1e-1;

KD = spdiags(K,0);
KD = spdiags(KD,0,n,n);

KB = K-KD;

resnorm = zeros(maxit,1);

figure(2);
U = zeros(n,1);
for ctr=1:maxit
%for ctr=maxit+1:2*maxit  
  U = KD\(f-KB*U);
  
  residual = K*U-f;
  
  if 0 % set to 1 if want to visualise  
    subplot(2,1,1);  
    plot(x,U,'b',x,u,'r')
    legend('Jacobi','slash');
    title(['iteration ',num2str(ctr)]);
    subplot(2,1,2)
    plot(x,residual);
    title('residual')
    drawnow
  end
  

  residual = norm(residual,2)./norm(f,2); % 2-norm of relative residual
  
  resnorm(ctr) = residual;
  
  if residual<toler
    break % leave the loop
  end
  
end


figure(3)
plot(resnorm,'.');


% analyse the properties of iteration matrix
T = KD\KB;
E = eig(full(T));
EE = eig(full(T'*T));
% spectral radius
disp('spectral radius of T:')
disp(max(abs(E)));
disp('spectral radius of T*T:')
disp(max(abs(EE)));










