% SpectralDiffFFTDemo.m illustrates periodic spectral differentiation with FFTs
% (cf. SpectralDiffDemo.m for the differentiation matrix implementation)
% for 2 different functions (one smooth, the other not)
% For complex y, delete "real" commands.
% Based on Trefethen: Spectral Methods in Matlab

clear all
clc
close all
set(gcf, 'units','normalized', 'pos',[.05 .25 .90 .65])

% Differentiation of a hat function:
N = 24; h = 2*pi/N; t = h*(1:N)';
y = max(0,1-abs(t-pi)/2);    % the hat function and
y_hat = fft(y);              % its discrete Fourier transform
yprime_hat = 1i*[0:N/2-1 0 -N/2+1:-1]' .* y_hat;   % differentiation in the Fourier domain
yprime = real(ifft(yprime_hat)); clf               % and inverse Fourier transform
subplot(2,2,1), plot(t,y,'.-','markersize',13)
axis([0 2*pi -.5 1.5]), grid on, title('function')
subplot(2,2,2), plot(t,yprime,'.-','markersize',13)
axis([0 2*pi -1 1]), grid on, title('spectral derivative')

% Differentiation of exp(sin(x)):
y = exp(sin(t)); yprime_exact = cos(t).*y;
y_hat = fft(y);
yprime_hat = 1i*[0:N/2-1 0 -N/2+1:-1]' .* y_hat;
yprime = real(ifft(yprime_hat));
subplot(2,2,3), plot(t,y,'.-','markersize',13)
axis([0 2*pi 0 3]), grid on
subplot(2,2,4), plot(t,yprime,'.-','markersize',13)
axis([0 2*pi -2 2]), grid on
error = norm(yprime-yprime_exact,inf);
text(2.2,1.4,['max error = ' num2str(error)])
