From: dirk@puk.ac.za (Dirk Laurie)
Newsgroups: sci.math.num-analysis
Subject: Re: Looking for numerical sum package
Date: 3 Oct 1995 06:26:29 GMT

M.C.W._van Rossum (vrossum@electron.phys.uva.nl) wrote:
: Hi, 

: I am looking for a package calculating sums like:
: sum_(i=1..infinity) f(i).
: ( f(i) of course decaing rapidly enough)

The best all-purpose algorithm is Levin's u transformation.  The Levin
transformations are special cases of the following:

Let the partial sums be s(n) and let a(n) be an auxiliary sequence tending
to zero.  Then 
s(n) = (s(n)/a(n)) / (1/a(n))
This is an infinite/infinite limit.  Regard the numerator and denominator
as functions of (1/n) and apply L'Hopital's rule.  You can't do it directly,
because you only know the functions at certain points, so use divided
differences instead of derivatives.

The u transformation is obtained using 
a(i) = i * f(i)

The following piece of Matlab 3.5 code is an implementation of the u
transformation with roundoff error estimate:

% LEVIN  y=levin(x,a) returns the diagonal of Levin's U-transformation on the 
%        (slowly convergent) sequence x, with shift a.
%        [y,e]=levin(x,a) also returns the estimated roundoff error in y
%        in terms of units in the last significant place.
%        levin(x) is the same as levin(x,0).
function [y,e]=levin(x,a)
n=length(x); if nargin<2, a=0; end
j=1:n; d=ones(j)./((j+a).*(x-[0,x(1:n-1)])); x=x.*d; t=abs(x);
for i=2:n, k=i:n; lk=(k+a).*(a+(1:n-i+1));
  x(k)=lk.*(x(k)-x(k-1)); t(k)=lk.*(t(k)+t(k-1)); d(k)=lk.*(d(k)-d(k-1)); 
end
y=x./d; e=t./abs(x);

The way to use it is: 
[y,e]=levin(cumsum(f))

The decision on which term of y to use (trade-off between accuracy of
approximation and build-up of roundoff error) is not so easy to automate.
I usually plot abs(diff(y)) and eps*e on the same graph and then it is
easy to make a decision.

Dirk Laurie