% Advection diffusion with small viscosity
%
%%%%%%%%%%%%%%%%%%%%%%%%%
%
% visc* u_xx - u_x = f
%
% u(0)=alpha, u(1)=beta, alpha=1, beta=3
%
% f=-1
%
% In this case (you can verify that) the exact solution is
%
% u=alpha+x+(beta-alpha-1)*(exp(x/visc)-1)./(exp(1/visc)-1);
% 
%%%%%%%%%%%%%%%%%%%%%%%%
%
% The solution of the problem is plotted against the numerical
% approximation (see below).
%
%%%%%%%%%%%%%%%%%%%%%%%%%
%
% You can run the program with different values of the viscosity parameter
% visc = 0.3, 0.1, 0.05
% and see that the "boundary layer" gets steaper and steaper. This fact
% creates difficulties for the numerical method.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%
% Try also to see what happens if you replace the backward difference
% operator D with central differences and modify the boundary term so that
% you get an overall second order method in h (the space discretization
% parameter). You'll see a bad oscillation near the boundary layer and only
% using h small enough (m bigger) this undesired effect will go away.
%
% This is related to what happens to u_t+au_x=0 (a=1): using central differences
% to approximate u_x and forward differences to approximate u_t you get a
% method that is always unstable (independently on h and k).
% But to fix it you can add a term corresponding to the central difference 
% discretization of u_xx ("add numerical viscosity") and get the Lax-Wendroff 
% method (see ch. 7.3).
%
% You can also try to solve
% u_t+au_x=visc* u_xx-f
% a=1, u(0)=alpha, u(1)=beta, u(x,0)=g(x)
% Using similar methods. You'll get a solution evolving from the initial
% function g(x) to the solution of visc*u_xx-au_x=f
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
T=1;
alpha=1;
beta=3;
%
m=100;
h=T/(m+1);
xx=[h:h:T-h]';
visc=0.005;
u=alpha+xx+(beta-alpha-1)*(exp(xx/visc)-1)./(exp(1/visc)-1);
plot([0;xx;T],[alpha;u;beta]);
e=ones(m,1);


A=1/(h^2)*spdiags([e -2*e e],-1:1,m,m);
%D=1/(2*h)*spdiags([-e e],[-1 1],m,m);
D=1/h*spdiags([-e e],-1:0,m,m);
F=-ones(m,1);
F(1,1)=F(1,1)-alpha*visc/(h^2)- alpha/(h);%alpha/(2*h);
F(end,1)=F(end,1)-beta*visc/(h^2);

Am=visc*A-D;
U=Am\F;
hold on
plot([0;xx;T],[alpha;U;beta],'r')
hold off