% SpectralAccuracyDemo.m illustrates the accuracy of periodic spectral differentiation
% for 4 different functions of increasing smoothness
% Based on Trefethen: Spectral Methods in Matlab

clear all
clc
close all
scrnsz = get(0, 'Screensize');
figure('Position', [5 10 scrnsz(3)-5 scrnsz(4)-60])

% Compute derivatives for various values of N:
Nmax = 50; E = zeros(3,Nmax/2-2);
for N = 6:2:Nmax;
    h = 2*pi/N; t = h*(1:N)';
    column = [0 .5*(-1).^(1:N-1).*cot((1:N-1)*h/2)]';
    D = toeplitz(column,column([1 N:-1:2]));
    y = abs(sin(t)).^3;                     % 3rd deriv in BV
    yprime = 3*sin(t).*cos(t).*abs(sin(t));
    E(1,N/2-2) = norm(D*y-yprime,inf);
    y = exp(-sin(t/2).^(-2));               % C-infinity
    yprime = .5*y.*sin(t)./sin(t/2).^4;
    E(2,N/2-2) = norm(D*y-yprime,inf);
    y = 1./(1+sin(t/2).^2);                 % analytic in a strip
    yprime = -sin(t/2).*cos(t/2).*y.^2;
    E(3,N/2-2) = norm(D*y-yprime,inf);
    y = sin(10*t); yprime = 10*cos(10*t);   % band-limited
    E(4,N/2-2) = norm(D*y-yprime,inf);
end

% Plot results:
titles = {'|sin(t)|^3','exp(-sin^{-2}(t/2))',...
       '1/(1+sin^2(t/2))','sin(10t)'}; clf
for iplot = 1:4 
    subplot(2,2,iplot)
    semilogy(6:2:Nmax,E(iplot,:),'.','markersize',12)
    line(6:2:Nmax,E(iplot,:))
    axis([0 Nmax 1e-16 1e3]), grid on
    set(gca,'xtick',0:10:Nmax,'ytick',(10).^(-15:5:0))
    xlabel N, ylabel error, title(titles(iplot))
end
