% Spectral differentiation compared
% with two finite difference schemes

clear;
clear all;
clear session;

k    = 0;
Nvec = 10:2:300;

for N=Nvec
    
    k = k + 1;
    h = 1.0 / N;
    f = 2.0 * pi * 1i * [0:N/2-1 0 -N/2+1:-1];
    
    x = linspace(0, 1 - 1/N, N);
    u = exp(-sin(pi*x).^(-2));
    v = 1 ./ (1 + sin(pi*x).^2);
    w = exp(sin(2.0*pi*x));
    
    du = pi * u .* sin(2.0*pi*x) ./ sin(pi*x).^4;
    dv = -2.0 * pi * sin(pi*x) .* cos(pi*x) .* v.^2;
    dw = 2.0 * pi * exp(sin(2.0*pi*x)) .* cos(2.0*pi*x);
    
    up1 = [u(2:N) u(1)];
    um1 = [u(N) u(1:N-1)];
    vp1 = [v(2:N) v(1)];
    vm1 = [v(N) v(1:N-1)];
    wp1 = [w(2:N) w(1)];
    wm1 = [w(N) w(1:N-1)];
    
    du1 = (up1 - u) / h;
    du2 = (up1 - um1) / (2*h);
    du3 = real(ifft(f .* fft(u)));
    dv1 = (vp1 - v) / h;
    dv2 = (vp1 - vm1) / (2*h);
    dv3 = real(ifft(f .* fft(v)));
    dw1 = (wp1 - w) / h;
    dw2 = (wp1 - wm1) / (2*h);
    dw3 = real(ifft(f .* fft(w)));
    
    eu1(k) = norm(du1-du,inf);
    eu2(k) = norm(du2-du,inf);
    eu3(k) = norm(du3-du,inf);
    ev1(k) = norm(dv1-dv,inf);
    ev2(k) = norm(dv2-dv,inf);
    ev3(k) = norm(dv3-dv,inf);
    ew1(k) = norm(dw1-dw,inf);
    ew2(k) = norm(dw2-dw,inf);
    ew3(k) = norm(dw3-dw,inf);
    
end;

fu = figure;
loglog(Nvec,eu1,'b',Nvec,eu2,'g',Nvec,eu3,'r');
legend('1st order','2nd order','spectral',3);
grid on;
axis([10 300 10^(-15) 10]);
xlabel('N');
ylabel('error');
title('exp(-sin(\pi x)^{-2})');

fv = figure;
loglog(Nvec,ev1,'b',Nvec,ev2,'g',Nvec,ev3,'r');
legend('1st order','2nd order','spectral',3);
grid on;
axis([10 300 10^(-15) 10]);
xlabel('N');
ylabel('error');
title('1 / (1 + sin(\pi x)^2)');

fw = figure;
loglog(Nvec,ew1,'b',Nvec,ew2,'g',Nvec,ew3,'r');
legend('1st order','2nd order','spectral',3);
grid on;
axis([10 300 10^(-15) 10]);
xlabel('N');
ylabel('error');
title('exp(sin(2 \pi x))');