%% Script for using Fourier spectral method to solve 
%
%       u" - u = 1/(1-0.6*cos(x)), 0 < x < L = 2*pi
% 
% Figure 1: 0<x<L, periodic BCs
% Figure 2: 0<x<L/2, homo Dirichlet BCs (using odd extension, 
%                                        induces discontinuity in f)
%
% Solution is plotted for N = 2^p, p = 1,..,5

  L = 2*pi;
  color = 'kbgmr';
  pmax = 5;
  pintp = 5;
  Nintp = 2^pintp;

%% Figure 1: 0<x<L, periodic BCs
  
  figure(1)
  clf
  delta = zeros(1,pmax);
  h = zeros(1,pmax);
  for p = 1:pmax
    N = 2^p;
    h(p) = L/N;
    M = [0:(N/2-1) (-N/2):(-1)];
    k = 2*pi*M/L;
    x = (0:(N-1))*h(p); % Fourier gridpoints
     
    f = 1./(1-0.6*cos(2*pi*x/L)); % RHS at gridpoints
           
    %    DFT, solve, and inverse DFT

    fhat = fft(f);
    u = real(ifft(-fhat./(k.^2+1)));
    
    uintp = interpft(u,Nintp);
    xintp = L*(0:Nintp - 1)/Nintp;
    subplot(2,3,2)
    plot(x/L,u,['x' color(p)],xintp/L,uintp,['-' color(p)])
    text(0.5,-1.3-0.1*p,['N=' int2str(N)],'Color',color(p))
    if (p==1)
      xlabel('x/L')
      ylabel('u')
      title('Fourier spectral u(x)')
      hold on
    else
      % Max-norm delta of h(p) - h(p-1) solution difference on common gridpoints     
      delta(p-1) = norm(u(1:2:N-1)-upm1,inf);
    end
    upm1 = u; % Define p-1 solution for next p
  end
  hold off

  % Plot f at the highest grid resolution used
  subplot(2,3,1)
  plot(x/L,f)
  xlabel('x/L')
  ylabel('f')
  title('f(x)')
 
  % Log-log plot of max-norm of h vs. h/2 solution difference
  subplot(2,3,3)
  loglog(h,delta,'x')
  xlabel('h')
  ylabel('\delta')
  title('Max-norm h/2-h soln diff')
  
 %% Test Dirichlet BCs u(0) = u(L/2) = 0
 
  figure(2)
  clf
  delta = zeros(1,pmax);
  h = zeros(1,pmax);
  for p = 1:pmax
    N = 2^p;
    h(p) = L/N;
    M = [0:(N/2-1) (-N/2):(-1)];
    k = 2*pi*M/L;
    x = (0:(N-1))*h(p); % Fourier gridpoints
    if (N>2)
      xh = x(2:N/2); 
      % RHS at interior gridpoints with 0 < x < L/2
      fh = 1./(1-0.6*cos(2*pi*xh/L)); 
      f  = [0 fh 0 -flip(fh)];  % Odd extension of f
    else
      f = nan(1,N); % If N/2=1, odd extn of f meaningless
    end
    
    %    DFT, solve, and inverse DFT

    fhat = fft(f);
    u = real(ifft(-fhat./(k.^2+1)));
    
    uintp = interpft(u,Nintp);
    xintp = L*(0:Nintp - 1)/Nintp;
    subplot(2,3,2)
    plot(x/L,u,['x' color(p)],xintp/L,uintp,['-' color(p)])
    text(0.7,-0.2-0.1*p,['N=' int2str(N)],'Color',color(p))
    if (p==1)
      xlabel('x/L')
      ylabel('u')
      title('Fourier spectral u(x)')
      hold on
    else
      % Max-norm delta of h(p) - h(p-1) solution difference on common gridpoints     
      delta(p-1) = norm(u(1:2:N-1)-upm1,inf);
    end
    upm1 = u; % Define p-1 solution for next p
  end
  hold off

  % Plot f at the highest grid resolution used
  subplot(2,3,1)
  plot(x/L,f)
  xlabel('x/L')
  ylabel('f')
  title('Odd extension of f_H(x)')
 
  % Log-log plot of max-norm of h vs. h/2 solution difference
  subplot(2,3,3)
  loglog(h,delta,'x')
  xlabel('h')
  ylabel('\delta')
  title('Max-norm h/2-h soln diff')
  hold on
  plot(h,0.1*h.^2,'r-')
  text(0.5*h(3),0.2*h(3).^2,'O(h^{-2})','Color','red')