% Spectral differentiation compared
% with two finite difference schemes

clear;
clear all;
clear session;

Nvec = 10:2:1000;
Mvec = 3:2:15;

for M=Mvec
    
    k = 0;
    
    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 = abs(sin(2.0*pi*x)).^M;
    
        du = 2.0 * pi * M * sin(2.0*pi*x) .* ...
             cos(2.0*pi*x) .* abs(sin(2.0*pi*x)).^(M-2);
    
        up1 = [u(2:N) u(1)];
        um1 = [u(N) u(1:N-1)];
        
        du1 = (up1 - u) / h;
        du2 = (up1 - um1) / (2*h);
        du3 = real(ifft(f .* fft(u)));
     
        eu1(k) = norm(du1-du,inf);
        eu2(k) = norm(du2-du,inf);
        eu3(k) = norm(du3-du,inf);
       
    end;

    figure;
    loglog(Nvec,eu1,'b',Nvec,eu2,'g',Nvec,eu3,'r');
    legend('1st order','2nd order','spectral',3);
    grid on;
    axis([10 1000 10^(-14) 10]);
    xlabel('N');
    ylabel('error');
    ts = ['|sin(2 \pi x)', sprintf('|^{%d}',M)];
    title(ts);
    
end;