% implicit.m: Solve the linear advection equation with 
%             an implicit scheme (BTCS or BTBS).
%
% Date:   January 25, 2011
% Author: John Stockie
%         Department of Mathematics
%         Simon Fraser University
%

function implicit( method, CFL, nx )

if nargin < 1,
  disp(  'Usage:  implicit( method, CFL, nx )' );
  disp( ['        where method is one of ''btcs'', ''btbs'''] );
  disp( ['Defaults: (''btcs'', 0.8, 100)'] );
  return;
end

if nargin < 2, CFL = 0.8; end;
if nargin < 3, nx  = 100; end;

a    = 1.0;   % advection speed
tend = 5;     % end time
xmax = 10;    % domain length [0,xmax]

dx = xmax/nx; % mesh spacing (constant)
x  = [0 : dx : xmax]';
ntprint = 50; % for printing
xex= [0 : xmax/1000 : xmax]';

% Calculate the time step (dt) needed to satisfy the CFL condition.
maxdt = CFL * dx / a;
nt    = ceil( tend / maxdt );
dt    = tend / nt;
sigma = a * dt / dx;

% Set up the initial solution values.
u0   = uexact(x, 0, a);
u    = u0;
unew = 0*u;

disp( ['    a = ', num2str(a)] );
disp( ['   dx = ', num2str(dx)] );
disp( ['sigma = ', num2str(sigma)] );

ntprint = min(nt, ntprint);
dtprint = tend / ntprint;

% Step the solution in time.
uall = zeros(nx+1,ntprint+1);
uall(:,1) = u0;
uerr0 = zeros(1,nt);
uerr1 = zeros(1,nt);
uerr2 = zeros(1,nt);

ee    = ones(nx+1,1);
if lower(method) == 'btcs', 
  A = spdiags([-sigma/2*ee, ee, sigma/2*ee], [-1:1], nx+1, nx+1);
else               % btbs
  A = spdiags([-sigma*ee, (1+sigma)*ee], [-1:0], nx+1, nx+1);
end
A(1,2)       = 0.0;
A(end,end-1) = 0.0;

ip = 1;
figure(1)

for i = 1 : nt,

  t = i*dt;
  uex = uexact(xex, t, a);
  
  unew = A \ u;
  % NOTE: no correction for boundary values here 
  % because BC is u(0,t) = 0.

  % Calculate the max-norm and L1-norm errors.
  uerr0(i) = max( abs(unew-uexact(x,t,a)) );
  uerr1(i) = sum( dx*abs(unew-uexact(x,t,a)) );
  uerr2(i) = sqrt( sum( dx*abs(unew-uexact(x,t,a)).^2 ) );
  
  % Plot the solution profiles.
  if t >= ip*dtprint,
    plot(x, unew)
    xlabel('x'), ylabel('u')
    title( ['Solution at time t=', num2str(t,'%9.4f')] )
    hold on, plot(xex,uex,'r:'), hold off
    grid on, shg
    pause(0.1)
    ip = ip + 1;
    uall(:,ip) = unew;
  end
  
  u = unew;

end

disp( ['U_max = ', num2str(max(u))] );
disp( ['L-1   = ', num2str(uerr1(end))] );
disp( ['L-2   = ', num2str(uerr2(end))] );
disp( ['L-inf = ', num2str(uerr0(end))] );

% Plot a sequence of profiles throughout the run.
figure(2)
nskip = 3;
plot(x,uall(:,1:nskip:end));
xlabel('x'), ylabel('u')
grid on, shg

% Plot the variation in error over time.
figure(3)
loglog(dt*[1:nt],[uerr0;uerr2;uerr1])
xlabel('t'), ylabel('Error')
legend('max-norm', 'L2-norm', 'L1-norm',0) 
axis tight
grid on, shg

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function uex = uexact( x, t, a )

global ictype;

u0shift = 2;
uex = 0.5*(1 + sin(2*pi*(x-1/4-u0shift-a*t))) .* ...
      (abs(x-0.5-u0shift-a*t)<0.5);
