% Demo quadrature of oscillatory function cos(1/x) on (0,1) via complex plane.
% Needs: gauss.m.  Barnett 5/22/14. Influences: "SIAM 100 Digit Challenge" book.
clear

f = @(z) exp(-1i./z);   % integrand. real part is what we want: cos(1/x)
figure; x = 0:1e-4:1; plot(x,real(f(x)),'-'); title('integrand Re f on (0,1)')

figure; x = -.5:0.002:1; [xx yy] = meshgrid(x); ff=f(xx+1i*yy);
imagesc(x,x,real(ff)); caxis([-2 2]); axis equal xy; title('f in complex plane')
hold on; plot([0 1], [0 0], '--.');  % desired integration path (osc)

z = @(x) x + 1i*(x-x.^2);        % change variable from z to x
x = 0:1e-3:1; plot(z(x),'--')    % add new integration path to plot

g = @(x) exp(-1i./(x + 1i*(x-x.^2))).*(1+1i-2i*x);  % g(x) dx = f(z(x)) dz(x)

figure; x = 0:1e-3:1; plot(x,real(g(x)),'-'); title('new integrand g (not osc)')

for N=30:30:180  % check convergence of usual quadr scheme for new integrand
  [x w] = gauss(N); x = x/2+1/2; w = w/2; disp(real(sum(g(x).*w'))) % 15 digits
end % (note Clenshaw-Curtis fails due to endpoint being a NaN)
