Singular Fourier-Pade approximation

Partial sums of Fourier terms for a function with jumps in value or derivative converge poorly, because of the Gibbs phenomenon. This file uses the Fourier coefficients, and locations of the singularities, to construct a different approximation that converges spectrally. For details, see T. A. Driscoll and B. Fornberg, Numerical Algorithms 26 (2001), pp. 77-92.

Example for f(x)=|x|, using 7 Fourier coefficients:

c = [pi/4 zeros(1,11)];
c(2:2:12) = -(2/pi)*(1:2:11).^(-2);
z0 = exp(1i*[-pi 0]);
[p,q,r] = padelog(c,z0);

% Make a plot:
x = linspace(-pi+10*eps,pi-10*eps,200); z = exp(1i*x);
pz = polyval(p(end:-1:1),z);
qz = polyval(q(end:-1:1),z);
rz{1} = polyval(r{1}(end:-1:1),z);
rz{2} = polyval(r{2}(end:-1:1),z);
fplus = ( pz + rz{1}.*log(1-z/z0(1)) + rz{2}.*log(1-z/z0(2)) ) ./ qz;
plot(x,abs(x),x,2*real(fplus),'k.')