altsum(a,method)

function s = altsum(a,method)
%ALTSUM Convergence acceleration of alternating series.
%	ALTSUM(A) where A is a vector, estimates the sum of the infinite 
%   alternating series A(1) - A(2) + A(3) - A(4) + ...
%
%   ALTSUM(A,METHOD) where A is a vector of length N and METHOD is a 
%   string selects one of the following linear acceleration methods:
%   'euler', Euler's method with error ~ 2^-N
%   'binom', a binomial weighting with error ~ 3^-N
%   'cheby', a weighting based on Chebyshev polynomials, error ~ 5.828^-N
%   Default method is 'cheby'.
%
%   EXAMPLES
%
%	% log(2) = 1 - 1/2 + 1/3 - 1/4 + ...
%	altsum(1./(1:20))
%
%	% pi = 4 - 4/3 + 4/5 - 4/7 + ...
%	altsum(4./(1:2:39))
%
%   % Euler-Mascheroni constant
%   k = 2:20;
%   altsum(log(k)./k)/log(2) + log(2)/2
%
%   % Divergent series
%   altsum(1:9,'euler')
%   altsum((1:9).^2,'euler')

%   Author: Jonas Lundgren <[e-mail]> 2009

%   2009-09-07  First published.
%   2009-09-16  Two acceleration methods added.

if nargin < 1, help altsum, return, end
if nargin < 2, method = 'cheby'; end

% Index vector
n = numel(a);
k = (1:n)';

% Select method
switch method
    case 'euler'
        % Euler's method ~ 2^-n
        p = 1/2;
        Q = (n-k+1)./k;
    case 'binom'
        % Binomial weighting ~ 3^-n
        p = 2/3;
        Q = (n-k+1)./k/2;
    otherwise
        % Chebyshev weighting ~ 5.828^-n
        p = 4/(3+sqrt(8));
        Q = (n-k+1)./k.*(n-k+1/2)./(2*n-k);
end

% Distribution
d = exp(cumsum([n*log(p); log(Q)]));        % d = cumprod([p^n; Q]);
w = cumsum(d);

% Alternate signs
a = a(:);
a(2:2:n) = -a(2:2:n);

% Reverse order
a = a(n:-1:1);

% Weighted sum
s = sum(w(1:n).*a)/w(n+1);