function mit18086_error_analysis
%MIT18086_ERROR_ANALYSIS
%    Applies the upwind, Lax-Friedrichs, and Lax-Wendroff
%    method to the transport equation u_t=u_x for a
%    sequence of step sizes for final time and Courant
%    number fixed. The norm of the errors is plotted.

% 02/2007 by Benjamin Seibold
% Feel free to modify for teaching and learning.
% ------------------------------------------------------------------
r = .9;                                       % fixed Courant number
n = round(10.^linspace(1,3,12)); % sequence of numbers of gridpoints
tf = .8;                                                % final time
type_ic = 2;                             % type of initial condition
% ------------------------------------------------------------------
clf
E = zeros(length(n),3);
for k = 1:length(n)
   x = linspace(-1,1,n(k))';
   h = x(2)-x(1);
   dt = r*h;
   u0 = f(x,type_ic);
   uf = f(x+ceil(tf/dt)*dt,type_ic);
   for i = 1:3, u{i} = u0; end
   I = eye(n(k));
   R = diag(ones(1,n(k)-1),1);
   name{1} = 'upwind';
   A{1} = (R-I)/h;
   name{2} = 'Lax-Friedrichs';
   A{2} = (R-R')/(2*h)+(h^2/dt/2)*(R-2*I+R')/h^2;
   name{3} = 'Lax-Wendroff';
   A{3} = (R-R')/(2*h)+(dt/2)*(R-2*I+R')/h^2;
   for i = 1:3
      M{i} = I+dt*A{i};
      u{i} = M{i}^ceil(tf/dt)*u{i};
      E(k,i) = norm(u{i}-uf,2)/sqrt(n(k));
   end
   subplot(1,2,1)
   plot(x,u0,'k:',x,uf,'k-',x,u{1},'b.-',x,u{2},'r.-',x,u{3},'m.-')
   axis([-1 1 -.1 1.1])
   title(sprintf('solution after %d time steps',ceil(tf/dt)))
   subplot(1,2,2)
   loglog(n(1:k),E(1:k,1),'b.-',n(1:k),E(1:k,2),'r.-',...
          n(1:k),E(1:k,3),'m.-')
   axis([min(n) max(n) min(min(E)) max(max(E))])
   grid on
   legend(name{1},name{2},name{3},1)
   title('error')
   pause(.1)
end
hold on
for i = 1:3
   loglog([.1 1]*n(end),[10^(1+(i==3)) 1]*E(end,i),'k--')
end
hold off

function y = f(x,type)
% initial condition function
if nargin<1, type = 1; end
switch type
case 1, y = exp(-40*(x-.5).^2);
case 2, y = max((1-8*(x-.5).^2),0).^2;
case 3, y = x>.2&x<.8;
end