% Numerically calculate soln. to advection eqn. dq/dt + dq/dx = 0 on
% domain 0<x<1 with periodic BCs and smooth ICs. Compare spectral method with
% RK4 time-differencing scheme and Courant numbers 0.25 and 0.5 with 
% Lax-Wendroff (LW) at various grid spacings. Plot 2-norm error comparison.

  L = 1;  % Domain size
  tf = 1; % Final time

  np = 4; % Number of dx's to try
  err_sp = zeros(2,np); % Initialize error matrix for spectral method
                        % 2 rows since we will consider 2 Courant numbers
  h_sp = zeros(1,np); % Initialize vector of grid spacings
  for p = 1:np
    N = 2^(p+2); % Number of gridpoints
    x = L*(0:(N-1))/N;  %  x-gridpoints [1xN]
    dx = L/N;
    M = [0:(N/2-1) (-N/2):(-1)];  
    k = 2*pi*M/L;  % Wavenumbers [1xN].  

    q0 = 1./(4+3*cos(2*pi*x/L));  % Initial condition
%    q0 = sign(0.25-abs(x-0.4));   % Initial condition
    for inu = 1:2
      nu = 0.25*inu;
      dt = nu*dx; % Timestep
      nt = round(tf/dt); % Number of timesteps to take

      qhat = fft(q0);  % Initial Fourier expansion coeffs [1xN].
      for it = 1:nt

        % March forward dqhat/dt = -Shat using RK4
        % where Shat is DFT of S(q) = dq/dx

        d1 = -dt*1i*k.*qhat; 
        d2 = -dt*1i*k.*(qhat + 0.5*d1);
        d3 = -dt*1i*k.*(qhat + 0.5*d2);
        d4 = -dt*1i*k.*(qhat + d3);
        qhat = qhat + (d1 + 2*d2 + 2*d3 + d4)/6; % New qhat
      end  % it
      q = ifft(qhat); % Numerical q at tf [1xN]
      err_sp(inu,p) = sqrt(dx)*norm(q-q0); % L^2 grid norm of soln. error at t=1
    end  % inu
    h_sp(p) = dx; % Grid spacing
  end

  % Lax-Wendroff with Courant number nu = 0.8

  nu = 0.8; 
  np = 4; % Number of dx's to try
  err_LW = zeros(1,np); % Initialize error matrix for spectral method
                        % 2 rows since we will consider 2 Courant numbers
  h_LW = zeros(1,np); % Initialize vector of grid spacings
  for p = 1:np
    N = 2^(p+2);
    x = L*(0:(N+1))/N;  %  x-gridpoints [1x(N+2)] - includes ghost cells.
    dx = L/N;
    dt = nu*dx; % Timestep
    nt = round(tf/dt); % Number of timesteps to take
    q0 = 1./(4+3*cos(2*pi*x/L));  % Initial condition
%    q0 = sign(0.25-abs(x-0.4));   % Initial condition
    q = q0;
    for it = 1:nt
      F(2:(N+2)) = q(1:(N+1)) + 0.5*(1 - nu)*(q(2:(N+2))-q(1:(N+1))); %Flux
      q(2:(N+1)) = q(2:(N+1)) - nu*(F(3:(N+2))-F(2:(N+1))); % Update interior
      q(1) = q(N+1); % Update boundary cells using periodic BCs
      q(N+2) = q(2);
    end
    err_LW(p) = sqrt(dx)*norm(q-q0);  
    h_LW(p) = dx;  
  end

  subplot(2,2,1)
  loglog(h_sp,err_sp(1,:),'bx',h_sp,err_sp(2,:),'mx',h_LW,err_LW,'rx')
  xlabel('dx')
  ylabel('2-norm of error')
  legend('Spectral, \nu = 0.25', 'Spectral, \nu = 0.5', 'LW, \nu = 0.8',4)
  title('Spectral/RK4 vs. LW method on q_t+q_x=0') 
  hold on

  % Add dx^2 and dx^4 fit lines

  c_sp = 20;
  plot(h_sp,c_sp*h_sp.^4,'b--')
  text(.07,.008,'dx^4','Color','b')
  c_LW = 20;   
  plot(h_LW,c_LW*h_LW.^2,'r--')
  text(.07,.3,'dx^2','Color','r')

  text(.1,1.5e-4,'q_0=1/(4+3cos(\pix))')
%  text(.1,1.5e-4,'q_0=square wave')
  hold off

  
