% p21.m - convergence of the Fourier method for linear elliptic equations
a=1;
% algebraic convergence for f1
Nvec = 2.^(3:12); 
clf, subplot(1,2,1);
for N=Nvec
    h=2*pi/N;  x=0:h:2*pi-h;
    ii=1:N;    kk=[0:N/2-1 N/2 -N/2+1:-1];  %kk=N/2-abs(ii-N/2-1);
    ue=sin(abs(x-pi)).^3;
    f=10*sin(abs(x-pi)).^3 - 6*sin(abs(x-pi));

    % three steps to solve the linear system	    
    v=fft(f, N);
    v=v./(a+kk.^2);
    u=ifft(v,N);

    % plot the L2 and H1 norm error
    uerr=fft(u-ue,N)/N;
    err0 = norm(abs(uerr), 2);                 % L^2 norm
    err1 = norm(abs(uerr).*sqrt(a+kk.^2), 2);  % H^1 norm
    loglog(N, err0, 'r+', 'markersize', 14), hold on
    loglog(N, err1, 'go', 'markersize', 14), hold on
end
grid on, xlabel N, ylabel error
title('The L^2(+) and H^1(o) convergence for f^1');
loglog(Nvec, Nvec.^(-2), '--');
text(100, 1e-5, 'N^{-2}', 'fontsize', 18);

% spectral convergence for f2
Nvec = 4:2:28; 
subplot(1,2,2);
for N=Nvec
    h=2*pi/N;  x=0:h:2*pi-h;
    ii=1:N;    kk=[0:N/2-1 N/2 -N/2+1:-1]; %kk=N/2-abs(ii-N/2-1);
    ue=exp(sin(x));
    f=exp(sin(x)).*(sin(x).^2+sin(x));

    % three steps to solve the linear system	    
    v=fft(f, N);
    v=v./(a+kk.^2);
    u=ifft(v,N);

    % plot the L2 and H1 norm error
    uerr=fft(u-ue,N)/N;
    err0 = norm(abs(uerr), 2);                 % L^2 norm
    err1 = norm(abs(uerr).*sqrt(a+kk.^2), 2);  % H^1 norm
    semilogy(N, err0, 'r+', 'markersize', 14), hold on
    semilogy(N, err1, 'go', 'markersize', 14), hold on
end
grid on, xlabel N, ylabel error
title('The L^2(+) and H^1(o) convergence for f^2');
semilogy(Nvec, exp(-Nvec*1.5), '--');
text(12, 1e-12, 'exp(-1.5N)', 'fontsize', 16);

print(gcf, '-depsc', 'p21.eps');
