% Solution and error of the transport equation
%
%    u_t + a*u_x = 0
%
% x \in I=[a,b]  and  t \in [0,T] with the 
% initial data  u(x,0) = u_0(x) with different methods
%   Euler forward for time / cebtral for space
%   Lax-Friedrichs
%   Lax-Wendroff
%   Upwind
% These methods are presented for example in 
% Quarteroni/Sacco/Saleri: Numerische Mathematik 2

% Author     : Stefan Hüeber
% Date       : January 30, 2003
% Institution: University of Stuttgart,
%              Institut for Applied Analysis and Numerical Mathematics,
%              Numerical Mathematics for super computers  
% Version    : 1.0
% Minor modifications by Greg Fasshauer (Illinois Inst. of Technology) 05/02/03


clear all; close all;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Parameters 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Initial data
%u0 = inline('sin(pi*x)');
u0 = inline('atan((0.5-x)/0.002)');
% exact solution for computing the error
%u0T = 'sin(pi*(x-a*T))';
u0T = 'atan((0.5-x+a*T)/0.002)';

% velocity
a = 2;

% Intervall 
T = 1;
I = [0 1];

% Number of time steps
TT(1) = 700;
% Number of space steps
XX(1) = 100;


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% menu
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
flag = menu('Please select a method',...
						'Euler-forward for time / central for space',...
						'Lax-Friedrichs',...
						'Lax-Wendroff',...
						'Upwind',...
						'Errorplot');

disp('Please wait ...');


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
switch flag
 case 1
	errsteps = 1;
 case 2
	errsteps = 1;
 case 3
	errsteps = 1;
 case 4
	errsteps = 1;
 case 5
	% Discretization for first step of the error
	TT(1) = 5;
	XX(1) = 2;
	% Number of refinement steps for error plot
	errstepseuler = 5;
	errsteps = 7;
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Euler forward for time / central for space
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
TTeuler(1) = TT(1);
XXeuler(1) = XX(1);
if flag==1
	errstepseuler = 1;
end
if  flag==1 | flag==5
	for k = 1:errstepseuler
		NT = TTeuler(k);
		NX = XXeuler(k);
		
		dt = T/NT;
		dx = (I(2)-I(1))/NX;
		lambda = dt/dx;
		
		% Initial conditions
		SOL = zeros(NX+2*NT+1,2);
		SOLPLOT = zeros(NX+1,NT+1);
		for j = 1:NX+2*NT+1
			SOL(j,1) = u0(I(2)+(NT+1-j)*dx);
		end
		SOLPLOT(:,1) = SOL(NT+1:NT+NX+1,1);
		
		% Finite Difference-Scheme: Euler-forward in time / central in space
		for n = 1:NT
			for j = 1+n:NX+2*NT+1-n
				SOL(j,2) = SOL(j,1) - 0.5*lambda*a*(SOL(j-1,1)-SOL(j+1,1));
			end
		SOL(:,1) = SOL(:,2);
		SOLPLOT(:,n+1) = SOL(NT+1:NT+NX+1,1);
		end
	
		% Computation of the L2-error for t=T
		xx = linspace(I(1),I(2),NX+1);
		[L2errEuler(k)] = L2error(xx,dx,SOL(NT+1:NT+NX+1,1),u0T,'0*x',a,T);
		
		if flag~=5
			% Visualization
			time = ones(NX+1,1)*linspace(0,T,NT+1);
			space = linspace(I(2),I(1),NX+1)'*ones(1,NT+1);
			%surf(time,space,SOLPLOT,'EdgeColor','none');
			%xlabel('Zeit t','FontSize',16);
			%ylabel('Ort x','FontSize',16);
			surf(space,time,SOLPLOT,'EdgeColor','none');
			xlabel('Space x','FontSize',16);
			ylabel('Time t','FontSize',16);
			pause(0.5);
		end
		
		if k~=errstepseuler
			XXeuler(k+1) = 2*XXeuler(k);
			TTeuler(k+1) = 4*TTeuler(k);
		end
	end
end
	



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Lax-Friedrichs
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
TTlf(1) = TT(1);
XXlf(1) = XX(1);
if flag==2 | flag==5
	for k = 1:errsteps
		NT = TTlf(k);
		NX = XXlf(k);
		
		dt = T/NT;
		dx = (I(2)-I(1))/NX;
		lambda = dt/dx;
		
		% Initial conditions
		SOL = zeros(NX+2*NT+1,2);
		SOLPLOT = zeros(NX+1,NT+1);
		for j = 1:NX+2*NT+1
			SOL(j,1) = u0(I(2)+(NT+1-j)*dx);
		end
		SOLPLOT(:,1) = SOL(NT+1:NT+NX+1,1);
		
		% Finite Difference-Scheme: Lax-Friedrichs
		for n = 1:NT
			for j = 1+n:NX+2*NT+1-n
				SOL(j,2) = 0.5*(SOL(j-1,1)+SOL(j+1,1)) - 0.5*lambda*a*(SOL(j-1,1)-SOL(j+1,1));
			end
			SOL(:,1) = SOL(:,2);
		SOLPLOT(:,n+1) = SOL(NT+1:NT+NX+1,1);
		end
		
		% Computation of the L2-error for t=T
		xx = linspace(I(1),I(2),NX+1);
		[L2errlf(k)] = L2error(xx,dx,SOL(NT+1:NT+NX+1,1),u0T,'0*x',a,T);
		
		if flag~=5
			% Visualization
			time = ones(NX+1,1)*linspace(0,T,NT+1);
			space = linspace(I(2),I(1),NX+1)'*ones(1,NT+1);
			%surf(time,space,SOLPLOT,'EdgeColor','none');
			%xlabel('Zeit t','FontSize',16);
			%ylabel('Ort x','FontSize',16);
			surf(space,time,SOLPLOT,'EdgeColor','none');
			xlabel('Space x','FontSize',16);
			ylabel('Time t','FontSize',16);
			pause(0.5);
		end
			
		if k~=errsteps
			XXlf(k+1) = 2*XXlf(k);
			TTlf(k+1) = 2*TTlf(k);
		end
	end
end



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Lax-Wendroff
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
TTlw(1) = TT(1);
XXlw(1) = XX(1);
if flag==3 | flag==5
	for k = 1:errsteps
		NT = TTlw(k);
		NX = XXlw(k);
		
		dt = T/NT;
		dx = (I(2)-I(1))/NX;
		lambda = dt/dx;
		
		% Initial conditions
		SOL = zeros(NX+2*NT+1,2);
		SOLPLOT = zeros(NX+1,NT+1);
		for j = 1:NX+2*NT+1
			SOL(j,1) = u0(I(2)+(NT+1-j)*dx);
		end
		SOLPLOT(:,1) = SOL(NT+1:NT+NX+1,1);
		
		% Finite Difference-Scheme: Lax-Wendroff
		for n = 1:NT
			for j = 1+n:NX+2*NT+1-n
				SOL(j,2) = SOL(j,1) - 0.5*lambda*a*(SOL(j-1,1) - SOL(j+1,1))...
						+ 0.5*(lambda*a)^2*(SOL(j-1,1)-2*SOL(j,1)+SOL(j+1,1));
			end
			SOL(:,1) = SOL(:,2);
		SOLPLOT(:,n+1) = SOL(NT+1:NT+NX+1,1);
		end
		
		% Computation of the L2-error for t=T
		xx = linspace(I(1),I(2),NX+1);
		[L2errlw(k)] = L2error(xx,dx,SOL(NT+1:NT+NX+1,1),u0T,'0*x',a,T);
		
		if flag~=5
			% Visualization
			time = ones(NX+1,1)*linspace(0,T,NT+1);
			space = linspace(I(2),I(1),NX+1)'*ones(1,NT+1);
			%surf(time,space,SOLPLOT,'EdgeColor','none');
			%xlabel('Zeit t','FontSize',16);
			%ylabel('Ort x','FontSize',16);
			surf(space,time,SOLPLOT,'EdgeColor','none');
			xlabel('Space x','FontSize',16);
			ylabel('Time t','FontSize',16);
			pause(0.5);
		end
			
		if k~=errsteps
			XXlw(k+1) = 2*XXlw(k);
			TTlw(k+1) = 2*TTlw(k);
		end
	end
end



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Upwind
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
TTuw(1) = TT(1);
XXuw(1) = XX(1);
if flag==4 | flag==5
	for k = 1:errsteps
		NT = TTuw(k);
		NX = XXuw(k);
		
		dt = T/NT;
		dx = (I(2)-I(1))/NX;
		lambda = dt/dx;
		
		% Initial conditions
		SOL = zeros(NX+2*NT+1,2);
		SOLPLOT = zeros(NX+1,NT+1);
		for j = 1:NX+2*NT+1
			SOL(j,1) = u0(I(2)+(NT+1-j)*dx);
		end
		SOLPLOT(:,1) = SOL(NT+1:NT+NX+1,1);
		
		% Finite Difference-Scheme: Lax-Wendroff
		for n = 1:NT
			for j = 1+n:NX+2*NT+1-n
				SOL(j,2) = SOL(j,1) - 0.5*lambda*a*(SOL(j-1,1) - SOL(j+1,1))...
						+ 0.5*lambda*abs(a)*(SOL(j-1,1)-2*SOL(j,1)+SOL(j+1,1));
			end
			SOL(:,1) = SOL(:,2);
		SOLPLOT(:,n+1) = SOL(NT+1:NT+NX+1,1);
		end
		
		% Computation of the L2-error for t=T
		xx = linspace(I(1),I(2),NX+1);
		[L2errUW(k)] = L2error(xx,dx,SOL(NT+1:NT+NX+1,1),u0T,'0*x',a,T);
		
		if flag~=5
			% Visualization
			time = ones(NX+1,1)*linspace(0,T,NT+1);
			space = linspace(I(2),I(1),NX+1)'*ones(1,NT+1);
			%surf(time,space,SOLPLOT,'EdgeColor','none');
			%xlabel('Zeit t','FontSize',16);
			%ylabel('Ort x','FontSize',16);
			surf(space,time,SOLPLOT,'EdgeColor','none');
			xlabel('Space x','FontSize',16);
			ylabel('Time t','FontSize',16);
			pause(0.5);
		end
		
		if k~=errsteps
			XXuw(k+1) = 2*XXuw(k);
			TTuw(k+1) = 2*TTuw(k);
		end
	end
end

if flag == 5
	% Errorplot	
	loglog( XXeuler.*TTeuler, 20*((1./XXeuler).^2 + 1./ ...
																									 TTeuler),'r--','linewidth',3);
	hold on; box on;
	loglog( XXlf.*TTlf, 1.1*(((1./XXlf).^2)./(1./TTlf) + (1./TTlf) + (1./XXlf)),'g--','linewidth',3);
	loglog( XXlw.*TTlw, 3.2*((1./XXlw).^2 + (1./TTlw).^2),'b--', ...
					'linewidth',3);
	loglog( XXuw.*TTuw, 1.3*((1./XXuw) + (1./TTuw)),'m--', ...
					'linewidth',3);
	loglog(XXeuler.*TTeuler,L2errEuler,'ro-','Linewidth',3, ...
				 'MarkerFaceColor','r');
	
	loglog(XXlf.*TTlf,L2errlf,'go-','Linewidth',3, ...
				 'MarkerFaceColor','g');
	loglog(XXlw.*TTlw,L2errlw,'bo-','Linewidth',3, ...
				 'MarkerFaceColor','b');
	loglog(XXuw.*TTuw,L2errUW,'mo-','Linewidth',3, ...
				 'MarkerFaceColor','m');
	title('L2-error for t=T','FontSize',24,'Color','r');
	xlabel('N_t*N_x','FontSize',16,'Color','k');
	ylabel('Fehler','FontSize',16,'Color','k');
	set(gca,'FontSize',14);
	set(gca,'xlim',[min(XXlf.*TTlf), max(XXlf.*TTlf)],'ylim',[0 10]);
		legend('O(\Delta t + \Delta x^2)',...
					 'O(\Delta x^2 / \Delta t + \Delta t + \Delta x)',...
					 'O(\Delta t^2 + \Delta x^2)',...
					 'O(\Delta t + \Delta x)',...
					 'Euler / central','Lax-Friedrichs','Lax-Wendroff','Upwind',3);
end
