% Script FS_heat :  Fourier Spectral method on heat eqn with Dirichlet
% and Neumann BCs
%
%   q_t = q_xx   0 < x < 1
%   q(x,0) = exp(-(x - 0.5).^2/(2*sigma^2)), sigma = 0.1
%   BCs:   Dirichlet  q(1,t) = q(0,t) = 0
%       or Neumann    q_x(1,t) = q_x(0,t) = 0
%
% Use N modes on an extended periodic domain [0,2] with grid spacing dx =
% 2/N.   Output is plot of solution at times 0, T=0.1.
%=======================================================================
  
  % Input parameters

  L = 1;
  sigma = 0.1; % Initial Gaussian pulse width
  nu = 0.5;  % nu = dt/dx^2; a nondimensional measure of timestep dt
  
  % Loop over grid spacings using index p such that N = 2^p.  
  % For each N, plot the solution and the 2-norm 
  %    eps_(N/2) = ||y_N(x) - y_(N/2)(x)|| 
  % computed using (1) the shared gridpoints and (2)
  % the shared Fourier modes.
  
  pmin = 3;
  pmax = 6;
  Nmin = 2^pmin;
  T = 2*nu*(2*L/Nmin)^2;  % Final time (= 4 timesteps for smallest N)
  color = 'rmbg';
  sym = 'o*+x';
  
  for p = pmin:pmax
    if(p>pmin)
      udold = ud;
      udhatold = udhat;
      unold = un;
      unhatold = unhat;
    end
   
    nx = 2^p;
    
    % Derived parameters
  
    dx = 2*L/nx;
    dt = nu*dx^2;
    nt = round(T/dt);
    
    % Initial condition
    
    x = dx*(0:(nx/2));
    u0 = exp(-(x - 0.5).^2/(2*sigma^2));
    ud0 = [0 u0(2:nx/2) 0 -u0(nx/2:-1:2)]; % Odd periodic extn
    un0 = [u0(1:(nx/2+1))  u0(nx/2:-1:2)]; % Even periodic extn
    xx = dx*(0:(nx-1));
    k = [0:(nx/2-1) -nx/2:-1]*2*pi/(2*L);

    udhat = fft(ud0);
    unhat = fft(un0);
    for it = 1:nt
      udhat = udhat.*(1 - 0.5*dt*k.^2)./(1 + 0.5*dt*k.^2);
      unhat = unhat.*(1 - 0.5*dt*k.^2)./(1 + 0.5*dt*k.^2);
    end
    ud = ifft(udhat);
    un = ifft(unhat);

    % Plot solutions on periodically extended grid
    
    subplot(2,2,1)
    plot(xx,ud,[color(p-pmin+1) sym(p-pmin+1)])
    hold on
    subplot(2,2,2)
    plot(xx,un,[color(p-pmin+1) sym(p-pmin+1)])
    hold on

    if (p>pmin)

      subplot(2,2,3)

      % Difference of ud_N and ud_N/2 at their common gridpoints xj/dx = 0,2,4
       
      errd = ud(1:2:(nx-1))-udold(1:nx/2);
      errdrms = norm(errd)/sqrt(nx/2);  % rms error (= grid-norm error/sqrt(L))
 
      % Difference of uhat_N/nx and udhat_N/2/(nx/2).  By Parseval's theorem,
      % 2-norm of this should nearly match the gridpoint-derived rms error.
       
      errdhat = udhat/nx ...
        - [udhatold(1:nx/4) zeros(1,nx/2) udhatold((nx/4+1):nx/2)]/(nx/2);
      errdhatrms = norm(errdhat);
      loglog(nx/2,errdrms,[color(p-pmin+1) 'o']) 
      hold on
      loglog(nx/2,errdhatrms,[color(p-pmin+1) 'x']) 

      subplot(2,2,4)

      % Difference of un_N and un_N/2 at their common gridpoints xj/dx = 0,2,4
       
      errn = un(1:2:(nx-1))-unold(1:nx/2);
      errnrms = norm(errn)/sqrt(nx/2);  % rms error (= grid-norm error/sqrt(L))
 
      % Difference of uhat_N/nx and unhat_N/2/(nx/2).  By Parseval's theorem,
      % 2-norm of this should nearly match the gridpoint-derived rms error.
       
      errnhat = unhat/nx ...
        - [unhatold(1:nx/4) zeros(1,nx/2) unhatold((nx/4+1):nx/2)]/(nx/2);
      errnhatrms = norm(errnhat);
      loglog(nx/2,errnrms,[color(p-pmin+1) 'o']) 
      hold on
      loglog(nx/2,errnhatrms,[color(p-pmin+1) 'x']) 
    end

  end
  
  % Add u(x,0) to subplots of ud and un
    
  subplot(2,2,1)
  plot(xx,ud0,'k-')
  xlabel('x')
  ylabel('u')
  ylim([-1.1 1.1])
%  legend('N=8','N=16','N=32','N=64','t=0','Location','EastOutside')
  title('FS on u_t = u_{xx}, Dirichlet BCs')
  plot([1 1],[-1.1 1.1],'k--')
  hold off
  
  subplot(2,2,2)
  plot(xx,un0,'k-')
  xlabel('x')
  ylabel('u')
  ylim([-1.1 1.1])
  legend('N=8','N=16','N=32','N=64','t=0','Location','SouthEast')
  title('FS on u_t = u_{xx}, Neumann BCs')
  plot([1 1],[-1.1 1.1],'k--')
  hold off

  % Finish up error subplot, including adding a N^-4 best-fit line
  
  subplot(2,2,3)
  N = 2.^(pmin:(pmax-1));  
  Nm4fit = errdrms*(N/N(end)).^-4;
  loglog(N,Nm4fit,'k--')
  text(N(end),10*errdrms,'N^{-4} fit')
  xlabel('N')
  ylabel('\epsilon_N')
%  legend('Gridpoint','Fourier')
  title('rms |y_{2N}-y_N|, Dirichlet')
  hold off
  
  subplot(2,2,4)
  Nm4fit = errnrms*(N/N(end)).^-4;
  loglog(N,Nm4fit,'k--')
  text(N(end),10*errnrms,'N^{-4} fit')
  xlabel('N')
  ylabel('\epsilon_N')
  legend('Gridpoint','Fourier')
  title('rms |y_{2N}-y_N|, Neumann')
  hold off
