% Several methods for the advection equation u_t + a u_x = 0 are compared eps
% controls the multiple of the Laplacian that is added to the centered finite
% difference operator in space.  this introduces some numerical dissipation
% that smoothes out high frequencies. Different values of eps correspond to
% different methods (Forward Euler, Lax-Friedrichs, upwind and Lax-Wendroff).
% Of course it is more efficient to write the schemes directly without using
% matrices!  The leapfrog method is also implemented, however the first for
% loop has to be commented out and the last for loop uncommented.  See "Finite
% Difference Methods for Ordinary and Partial Differential Equations", R.
% LeVeque.
 

n=100;
x = linspace(0,1,n)';
utrue = @(x)  max(1-abs(6*(mod(x,1)-1/2)),0) + max(1-abs(20*(mod(x,1)-0.2)),0);
u = utrue(x);

a = 2;
e = ones(n,1); h = 1/n;
A = -(a/2/h)*spdiags([-e,e],[-1,1],n,n);
A(1,n)=a/2/h; A(n,1)=-a/2/h;

B = (1/h^2)*spdiags([e,-2*e,e],-1:1,n,n);
B(1,n)=1/h^2; B(n,1)=1/h^2;

k = 0.9*h/abs(a);

%eps = 0; % Euler (unstable)
%eps=h^2/2/k; % Lax-Friedrichs
eps=abs(a)*h/2;   % upwind
%eps=-a*h/2;   % upwind with wrong sign
%eps=a^2*k/2; % Lax-Wendroff
Aeps = A + eps*B;
figure(1); t=linspace(0,2*pi);
lambda=eig(full(k*Aeps));
plot(real(lambda),imag(lambda),'*',-1+cos(t),sin(t),'r-'); axis equal
fprintf('h=%f, k=%d, |a*k/h| = %f\n',h,k,abs(a*k/h));

figure(2)
% Forward Euler, Lax-Friedrichs, upwind and Lax-Wendroff
for i=1:n,
 u = u + k*Aeps*u;
 plot(x,u,x,utrue(x-a*k*i)); axis([0 1 -0.5 1.5]);
 title(sprintf('t=%d\n',k*i));
 pause(0.1);
end;

%% Leapfrog
%uol=u;
%u = uol + k*A*uol;
%for i=1:1000,
% unew = uol + 2*k*A*u;
% uol=u;
% u=unew;
% plot(x,u,x,utrue(x-a*k*i)); axis([0 1 -0.5 1.5]);
% title(sprintf('t=%d\n',k*(i-1)));
% pause(0.1);
%end;