% p11f.m - Comparison among the FDM, FEM, and Fourier spectral method
%          in solving the periodic heat equation
clear;
u0=inline('exp(sin(x))','x');

% Use the spectral procedure to calculate a solution as exact solution
function  [u, x] = EvalSolEx(uinit, N, t)
    h=2*pi/N;  x=(0:h:2*pi-h)';
    k=[0:N/2-1 N/2 -N/2+1:-1]'; 
    ui=feval(uinit, x);
    v=fft(ui, N); 
    v=v.*exp(-k.^2 * t);
    u=real(ifft(v,N)); 
end

function  errsp = spm(uinit, Nvec, t)
   errsp = []; 
   [ue, xe] = EvalSolEx(uinit, 2048, t);
   for N=Nvec
      [u, x] = EvalSolEx(uinit, N, t);
      solu = interp1(xe, ue, x, 'spline');
      err = norm(u-solu, 'inf');
      errsp = [errsp err]; 
   end 
end

% algebraic convergence for f1
function errfd = fdm(uinit, Nvec, dt, nt)
   errfd = [];
   [ue, xe] = EvalSolEx(uinit, 1024, dt*nt);
   for N=Nvec
       h =2*pi/N; x=(0:h:2*pi-h)'; a=dt/h^2;
       u0 = feval(uinit, x);
       D = toeplitz([2 -1 zeros(1,N-3) -1]);
       un = (eye(N)+a*D)^(-nt)*u0;
       solu = interp1(xe, ue, x, 'spline');
       err = norm(un-solu, 'inf');
       errfd = [errfd err]; 
   end
end

function errfe = fem(uinit, Nvec, dt, nt)
   errfe = [];
   [ue, xe] = EvalSolEx(uinit, 1024, dt*nt);
   for N=Nvec
       h = 2*pi/N; x=(0:h:2*pi-h)'; a=dt/h^2;
       u0 = feval(uinit, x);
       M = toeplitz([2/3 1/6 zeros(1,N-3) 1/6]);
       S = toeplitz([2 -1 zeros(1,N-3) -1]);
       un = (eye(N)+a*(M\S))^(-nt)*u0;
       solu = interp1(xe, ue, x, 'spline');
       err = norm(un-solu, 'inf');
       errfe = [errfe err]; 
   end
end

% algebraic convergence for FEM and FDM
figure('Position',[0 0 1000 500] ), subplot(1,2,1); 
Nvec = 2.^(3:9); dt = .5e-4; nt = 10000;
errfd = fdm(u0, Nvec, dt, nt);
errfe = fem(u0, Nvec, dt, nt);
loglog(Nvec, errfd, 'r+', 'markersize', 14), hold on
loglog(Nvec, errfe, 'bo', 'markersize', 14), hold on
title('The maximum error of FD(+) and FE(o) methods');
grid on, xlabel N, ylabel error
loglog(Nvec, Nvec.^(-2), '--');
text(100, 1e-5, 'N^{-2}', 'fontsize', 18);

% spectral convergence for spectral method
subplot(1,2,2);
Nvec = 4:2:16; 
errsp= spm(u0, Nvec, dt*nt);
semilogy(Nvec, errsp, 'r+', 'markersize', 14), hold on
grid on, xlabel N, ylabel error
title('The maximum error of the Fourier spectral method');
semilogy(Nvec, 1e3*exp(-Nvec*3), '--');
text(8, 1e-16, 'exp(-3N)', 'fontsize', 18);


