% Script to solve HW2 p8: Trapezoidal time and centered space
% differencing on advection equation du/t + c*du/dx = 0, u(0)=u(1),
% u(x,0) = sin(2*pi*x).

  % Parameters of problem
  c = 1;  % Advection speed
  L = 1;  % Domain width
  mu = 1; % Courant number
  T = 1;  % Final time
  
  % Loop over grid spacings
  
  prange = 1:4;
  np = length(prange);
  err = zeros(1,np);
  dx = 0.1*2.^-prange;
  
  for p = 1:np
    nx = round(L/dx(p)); % number of gridpoints 
    dt = mu*dx(p)/c;
    nt = round(T/dt);
    x = (1:nx)'*dx(p);   % x (as a column vector)
          
    % Set up numerical method
    % u(j,n+1) - u(j,n) = -0.25*(c*dt/dx)*(  u(j+1,n+1) + u(j+1,n)
    %                                      - u(j-1,n+1) - u(j-1,n))
    % in the form of an implicit matrix solve Mu(n+1) = rhs
    % The implicit solve makes it easier to treat periodic BCs
    % using indexing tricks rather than image points.
     
    jp1 = [2:nx 1];
    jm1 = [nx 1:(nx-1)];
    e = ones(nx,1);
    M = spdiags([-0.25*mu*e e 0.25*mu*e],-1:1,nx,nx);
    M(1,nx) = -0.25*mu;
    M(nx,1) = 0.25*mu;

    % Now timestep
    u = sin(2*pi*x);  % Initial u 
    for it = 1:nt
      rhs = u - 0.25*mu*(u(jp1)-u(jm1));
      u = M\rhs;
    end
    uT = sin(2*pi*x); % Exact soln u(x,T)
    err(p) = norm(u - uT)*sqrt(dx(p)); % Grid 2-norm of error
  end
  
  subplot(2,2,1)
  loglog(dx,err,'x',dx,err(p)*(dx/dx(np)).^2,'-')
  xlabel('\Delta x')
  ylabel('Grid 2-norm error at t=1')
  legend('Error','c\Delta x^2')
  title('Trapezoidal centered method')
  text(1.1e-3,0.5,'du/dt+du/dx=1')
  text(1.1e-3,0.25,'u(x,0)=sin2\pix')
  text(1.1e-3,0.125,'u(0,t)=u(1,t)')