%% Stationary Convection-Diffusion Equation
% This example computes a finite-difference-solution for a 1-D linear
% second-order boundary value problem with Dirichlet boundary conditions.
% CWB 2011-2017

%% The boundary value problem
% The dependent variable $c=c(x)$ models the concentration of a quantity
% along the x-Axis. The equation models convective transport with velocity
% $w$, balanced by diffusion with diffusivity $D$.
% Differential equation:
%
% $$D c'' -  w c' = 0  \mbox{ for }  c=c(x),\quad 0<x<L $$
%
% Boundary conditions
%
% $$ c = c_0 \mbox{ at } x=0, \quad c = c_L \mbox{ at } x=L $$
%% Parameters of the Problem
L = 1;                  % length of computational domain
Pe = 140;                % dimensionless Peclet-Number w L / D
c_0 = 1;                % Dirichlet BC on both sides
c_L = 0;

%% Discretization parameters
n = 200;                 % Nuber of interior gridpoints
x = linspace(0,L,n+2);
h = L/(n+1);            % grid size
xi = x(2:n+1);          % interor grid points
Pgr = Pe*h/L            % grid Peclet number

%% set up matrix equation:
% we try various ways to set up the tridiagonal matrix
%  'cd' ... central difference for convection term
%  'bd' ... backward difference
scheme='bd';
switch scheme
    case 'cd'
        % central difference stencil
        % (-1/h^2) * [-1-Pgr/2, 2, -1+Pgr/2]
        cm = (-1-Pgr/2)*ones(n,1);
        c0 = (2)*ones(n,1);
        cp = (-1+Pgr/2)*ones(n,1);
    case 'bd'
        % backward difference stencil  (-1/h^2) * [-1-Pgr, 2+Pgr, -1]
        cm = (-1-Pgr)*ones(n,1);
        c0 = (2+Pgr)*ones(n,1);
        cp = (-1)*ones(n,1);
end

A = spdiags([cm,c0,cp],[-1,0,1],n,n);

% right-hand side
b = zeros(n,1);
switch scheme
    case 'cd'
        b(1) = (1+Pgr/2)*c_0; % central
        b(n) = (1-Pgr/2)*c_L;
    case 'bd'
        b(1) = (1+Pgr)*c_0;    % backward
        b(n) = (1)*c_L;
end

ci = A\b;           % solution
c = [c_0;ci;c_L];   % solution and BC

%% Evaluation
figure(1)
xx=linspace(0,L,1000);
exact_c = c_0 + (exp(Pe*xx/L) - 1) / (exp(Pe) - 1) * (c_L-c_0);
plot(x,c,'-o',xx,exact_c)
title('Exact and Numerical Solution');

figure(2)
err = c' -(c_0 + (exp(Pe*x/L) - 1) / (exp(Pe) - 1) * (c_L-c_0));
plot(x,err')
title('Solution Error');

%% Abschneidefehler
% some calculations regarding truncation erron and stability bound of matrix.
% It turns out: these estimates are too pessimistic by some orders of magnitude! 
% trunc=1/12 *h^2*Pe^4*exp(Pe)/ (exp(Pe) - 1)
%
% ninv=norm(inv(A/h^2),inf)
% nA=norm(A/h^2,inf)
%
% trunc*ninv





