The functions described in this chapter accelerate the convergence of a series using the Levin u-transform. This method takes a small number of terms from the start of a series and uses a systematic approximation to compute an extrapolated value and an estimate of its error. The u-transform works for both convergent and divergent series, including asymptotic series.

These functions are declared in the header file gsl_sum.h.

Acceleration functions
======================

 -- Function: int gsl_sum_levin_u (const double * ARRAY, size_t ARRAY_SIZE, double * Q_NUM, double * Q_DEN, double * DQ_NUM, double * DQ_DEN, double * DSUM, double * SUM_ACCEL, double * SUM_PLAIN, double * PRECISION)

     This function takes the terms of a series in ARRAY of size ARRAY_SIZE and computes the extrapolated limit of the series using a Levin u-transform. The extrapolated sum is stored in SUM_ACCEL, with an estimate of the relative error stored in PRECISION. The actual term-by-term sum is returned in SUM_PLAIN. Additional working space must be provided in the arrays Q_NUM, Q_DEN, DSUM of size ARRAY_SIZE elements each and in DQ_NUM, DQ_DEN of size ARRAY_SIZE**2. The algorithm estimates both truncation error and round-off error to choose an optimal number of terms for the extrapolation.

Example of accelerating a series
================================

The following code calculates an estimate of $\zeta(2) = \pi^2 / 6$ \zeta(2) = \pi^2 / 6 using the series,

     \zeta(2) = 1 + 1/2^2 + 1/3^2 + 1/4^2 + ...

After N terms the error in the sum is O(1/N), making direct summation of the series converge slowly.

     #include <stdio.h>
     #include <gsl_math.h>
     #include <gsl_sum.h>

     #define N 20

     int
     main (void)
     {
       double t[N], qnum[N], qden[N], dsum[N], dqnum[N * N], dqden[N * N];
       double sum_accel, sum_plain, prec;
       double sum = 0;
       size_t n_used ;
       int n;

       const double zeta_2 = M_PI * M_PI / 6.0;

       /* terms for zeta(2) = \sum_{n=0}^{\infty} 1/n^2 */

       for (n = 0; n < N; n++)
         {
           double np1 = n + 1.0;
           t[n] = 1.0 / (np1 * np1);
           sum += t[n] ;
         }

       gsl_sum_levin_u_accel (t, N, qnum, qden, dqnum, dqden, dsum,
                              &sum_accel, &n_used, &sum_plain, &prec);

       printf("term-by-term sum = % .16f using %d terms\n", sum, N) ;
       printf("term-by-term sum = % .16f using %d terms\n", sum_plain, n_used) ;
       printf("exact value      = % .16f\n", zeta_2) ;
       printf("accelerated sum  = % .16f using %d terms\n", sum_accel, n_used) ;
       printf("estimated error  = % .16f\n", prec * sum_accel) ;
       printf("actual error     = % .16f\n", sum_accel - zeta_2) ;

       return 0;
     }

The ouput below shows that the Levin u-transform is able to obtain an estimate of the sum to 1 part in $10^10$ 10^10 using the first eleven terms of the series. The error estimate returned by the function is also accurate, giving the correct number of significant digits.

     bash$ ./a.out
     term-by-term sum =  1.5961632439130233 using 20 terms
     term-by-term sum =  1.5759958390005426 using 13 terms
     exact value      =  1.6449340668482264
     accelerated sum  =  1.6449340668166479 using 13 terms
     estimated error  =  0.0000000000508580
     actual error     = -0.0000000000315785

Note that a direct summation of this series would require $10^10$ 10^10 terms to achieve the same precision as the accelerated sum does in 13 terms.

References and Further Reading
==============================

The algorithms used by these functions are described in the following papers,

  T. Fessler, W.F. Ford, D.A. Smith, HURRY: An acceleration algorithm for scalar sequences and series ACM Transactions on Mathematical Software, 9(3):346--354, 1983. and Algorithm 602 9(3):355--357, 1983.

The theory of the u-transform was presented by Levin,

  D. Levin, Development of Non-Linear Transformations for Improving Convergence of Sequences, Intern. J. Computer Math. B3:371--388, 1973