%
% Finite difference discretization of the BVP
%
%     - diff  y''(x) + adv y'(x) + react y(x) = f(x)  x in (0,1),   
%                                        y(0) = g0,
%                                        y(1) = g1,
%
% Matthias Heinkenschloss
% Jan 14, 2020  (revised Jan 22, 2020)
%


set(0, 'defaultaxesfontsize',28,'defaultaxeslinewidth',3.0,...
       'defaultlinelinewidth',3.0,'defaultpatchlinewidth',3.0,...
       'defaulttextfontsize',28,'defaulttextinterpreter','latex');


%% problem data
diff  = 1.e-3;
adv   = 1;
react = 0;
g0    = 0;
g1    = 0;
  
  
f     = @(x) ones(size(x));
yex   = @(x) x - (exp((x-1)/diff)-exp(-1/diff))./(1 - exp(-1/diff));



%% finite difference discretization
upwind = 1; % =0 centered difference approx., =1 upwind

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

if( upwind == 0 )
    disp(' Using centered difference approximation')
    
    % system matrix
    % set subdiagonal
    c    = -(diff/h^2 + 0.5*adv/h)*ones(n,1);
    % set diagonal
    d    = (2*diff/h^2 + react)*ones(n,1);
    % set superdiagonal
    e    = -(diff/h^2 - 0.5*adv/h)*ones(n,1);
    % system matrix
    A = spdiags([c d e], -1:1, n, n);
    
    % right hand side
    b      = f(x(2:n+1));
    b(1)   = b(1) + (diff/h^2 + 0.5*adv/h)*g0;
    b(n)   = b(n) + (diff/h^2 - 0.5*adv/h)*g1;
    
    % FD  solution
    yfd = A\b;
    
    filename_sol = ['bvp_ex1_centered_n',int2str(n),'.eps'];
    
    
else
    disp(' Using upwind finite difference approximation')
    
    if( adv < 0 )
        error( 'adv < 0 not yet implemented')
    end
    % system matrix
    % set subdiagonal
    c    = -(diff/h^2 + adv/h)*ones(n,1);
    % set diagonal
    d    = (2*diff/h^2 + adv/h + react)*ones(n,1);
    % set superdiagonal
    e    = (-diff/h^2)*ones(n,1);
    % system matrix
    A = spdiags([c d e], -1:1, n, n);
    
    % right hand side
    b      = f(x(2:n+1));
    b(1)   = b(1) + (diff/h^2 + 0.5*adv/h)*g0;
    b(n)   = b(n) + (diff/h^2)*g1;
    
    % FD  solution
    yfd = A\b;
  
    filename_sol = ['bvp_ex1_upwind_n',int2str(n),'.eps'];

end

% plot solution
figure(1)
xx = (0:1/100:1);
plot(x,[g0;yfd;g1],'b-'); hold on
plot(xx,yex(xx),'r--'); hold off
xlabel('$x$')
ylabel('$y(x)$')
legend('FD approximation','exact solution','Location','NorthWest')
axis([-0.1,1.1,-0.1,1.1])
%print('-depsc',filename_sol)

%% discretization error analysis
for i = 1:14
    n     = 20*2^i - 1;
    h     = 1/(n+1);
    x     = (0:h:1)';
    hh(i) = h;

    % centered difference approximation
    % matrix
    c    = -(diff/h^2 + 0.5*adv/h)*ones(n,1);
    d    = (2*diff/h^2 + react)*ones(n,1);
    e    = -(diff/h^2 - 0.5*adv/h)*ones(n,1);
    A    = spdiags([c d e], -1:1, n, n);
    
    % right hand side
    b    = f(x(2:n+1));
    b(1) = b(1) + (diff/h^2 + 0.5*adv/h)*g0;
    b(n) = b(n) + (diff/h^2 - 0.5*adv/h)*g1;
    
    % FD  solution
    yfd = A\b;
    
    % error
    err1(i) = norm( yfd - yex(x(2:n+1)), 'inf');

    % upwind difference approximation

    % matrix
    c    = -(diff/h^2 + adv/h)*ones(n,1);
    d    = (2*diff/h^2 + adv/h + react)*ones(n,1);
    e    = (-diff/h^2)*ones(n,1);
    A    = spdiags([c d e], -1:1, n, n);
    
    % right hand side
    b    = f(x(2:n+1));
    b(1) = b(1) + (diff/h^2 + 0.5*adv/h)*g0;
    b(n) = b(n) + (diff/h^2)*g1;
    
    % FD  solution
    yfd = A\b;
    
    % error
    err2(i) = norm( yfd - yex(x(2:n+1)), 'inf');
end

figure(2)
loglog( hh, err1, 'r*-'); hold on
loglog( hh, err2, 'bo--'); hold on
loglog( [hh(9),hh(end)], 0.1*err1(9)/hh(9)^2*[hh(9)^2,hh(end)^2], 'k-'); hold on
text(hh(11),0.05*err1(11),'$h^2$'); hold on
loglog( [hh(9),hh(end)], 2*err2(9)/hh(9)*[hh(9),hh(end)], 'k-'); hold on
text(hh(11),8*err2(11),'$h$'); hold off
xlabel('$h$'); ylabel('error')
legend('centered','upwind','Location','SouthEast')
print('-depsc','bvp_err_ex1.eps')
