|
Packit |
67cb25 |
.. index::
|
|
Packit |
67cb25 |
single: acceleration of series
|
|
Packit |
67cb25 |
single: summation, acceleration
|
|
Packit |
67cb25 |
single: series, acceleration
|
|
Packit |
67cb25 |
single: u-transform for series
|
|
Packit |
67cb25 |
single: Levin u-transform
|
|
Packit |
67cb25 |
single: convergence, accelerating a series
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
*******************
|
|
Packit |
67cb25 |
Series Acceleration
|
|
Packit |
67cb25 |
*******************
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
The functions described in this chapter accelerate the convergence of a
|
|
Packit |
67cb25 |
series using the Levin :math:`u`-transform. This method takes a small number of
|
|
Packit |
67cb25 |
terms from the start of a series and uses a systematic approximation to
|
|
Packit |
67cb25 |
compute an extrapolated value and an estimate of its error. The
|
|
Packit |
67cb25 |
:math:`u`-transform works for both convergent and divergent series, including
|
|
Packit |
67cb25 |
asymptotic series.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
These functions are declared in the header file :file:`gsl_sum.h`.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
Acceleration functions
|
|
Packit |
67cb25 |
======================
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
The following functions compute the full Levin :math:`u`-transform of a series
|
|
Packit |
67cb25 |
with its error estimate. The error estimate is computed by propagating
|
|
Packit |
67cb25 |
rounding errors from each term through to the final extrapolation.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
These functions are intended for summing analytic series where each term
|
|
Packit |
67cb25 |
is known to high accuracy, and the rounding errors are assumed to
|
|
Packit |
67cb25 |
originate from finite precision. They are taken to be relative errors of
|
|
Packit |
67cb25 |
order :macro:`GSL_DBL_EPSILON` for each term.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
The calculation of the error in the extrapolated value is an
|
|
Packit |
67cb25 |
:math:`O(N^2)` process, which is expensive in time and memory. A faster
|
|
Packit |
67cb25 |
but less reliable method which estimates the error from the convergence
|
|
Packit |
67cb25 |
of the extrapolated value is described in the next section. For the
|
|
Packit |
67cb25 |
method described here a full table of intermediate values and
|
|
Packit |
67cb25 |
derivatives through to :math:`O(N)` must be computed and stored, but this
|
|
Packit |
67cb25 |
does give a reliable error estimate.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. type:: gsl_sum_levin_u_workspace
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
Workspace for a Leven :math:`u`-transform.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. function:: gsl_sum_levin_u_workspace * gsl_sum_levin_u_alloc (size_t n)
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
This function allocates a workspace for a Levin :math:`u`-transform of :data:`n`
|
|
Packit |
67cb25 |
terms. The size of the workspace is :math:`O(2n^2 + 3n)`.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. function:: void gsl_sum_levin_u_free (gsl_sum_levin_u_workspace * w)
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
This function frees the memory associated with the workspace :data:`w`.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. function:: int gsl_sum_levin_u_accel (const double * array, size_t array_size, gsl_sum_levin_u_workspace * w, double * sum_accel, double * abserr)
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
This function takes the terms of a series in :data:`array` of size
|
|
Packit |
67cb25 |
:data:`array_size` and computes the extrapolated limit of the series using
|
|
Packit |
67cb25 |
a Levin :math:`u`-transform. Additional working space must be provided in
|
|
Packit |
67cb25 |
:data:`w`. The extrapolated sum is stored in :data:`sum_accel`, with an
|
|
Packit |
67cb25 |
estimate of the absolute error stored in :data:`abserr`. The actual
|
|
Packit |
67cb25 |
term-by-term sum is returned in :code:`w->sum_plain`. The algorithm
|
|
Packit |
67cb25 |
calculates the truncation error (the difference between two successive
|
|
Packit |
67cb25 |
extrapolations) and round-off error (propagated from the individual
|
|
Packit |
67cb25 |
terms) to choose an optimal number of terms for the extrapolation.
|
|
Packit |
67cb25 |
All the terms of the series passed in through :data:`array` should be non-zero.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
Acceleration functions without error estimation
|
|
Packit |
67cb25 |
===============================================
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
The functions described in this section compute the Levin :math:`u`-transform of
|
|
Packit |
67cb25 |
series and attempt to estimate the error from the "truncation error" in
|
|
Packit |
67cb25 |
the extrapolation, the difference between the final two approximations.
|
|
Packit |
67cb25 |
Using this method avoids the need to compute an intermediate table of
|
|
Packit |
67cb25 |
derivatives because the error is estimated from the behavior of the
|
|
Packit |
67cb25 |
extrapolated value itself. Consequently this algorithm is an :math:`O(N)`
|
|
Packit |
67cb25 |
process and only requires :math:`O(N)` terms of storage. If the series
|
|
Packit |
67cb25 |
converges sufficiently fast then this procedure can be acceptable. It
|
|
Packit |
67cb25 |
is appropriate to use this method when there is a need to compute many
|
|
Packit |
67cb25 |
extrapolations of series with similar convergence properties at high-speed.
|
|
Packit |
67cb25 |
For example, when numerically integrating a function defined by a
|
|
Packit |
67cb25 |
parameterized series where the parameter varies only slightly. A
|
|
Packit |
67cb25 |
reliable error estimate should be computed first using the full
|
|
Packit |
67cb25 |
algorithm described above in order to verify the consistency of the
|
|
Packit |
67cb25 |
results.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. type:: gsl_sum_levin_utrunc_workspace
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
Workspace for a Levin :math:`u`-transform without error estimation
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. function:: gsl_sum_levin_utrunc_workspace * gsl_sum_levin_utrunc_alloc (size_t n)
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
This function allocates a workspace for a Levin :math:`u`-transform of :data:`n`
|
|
Packit |
67cb25 |
terms, without error estimation. The size of the workspace is
|
|
Packit |
67cb25 |
:math:`O(3n)`.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. function:: void gsl_sum_levin_utrunc_free (gsl_sum_levin_utrunc_workspace * w)
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
This function frees the memory associated with the workspace :data:`w`.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. function:: int gsl_sum_levin_utrunc_accel (const double * array, size_t array_size, gsl_sum_levin_utrunc_workspace * w, double * sum_accel, double * abserr_trunc)
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
This function takes the terms of a series in :data:`array` of size
|
|
Packit |
67cb25 |
:data:`array_size` and computes the extrapolated limit of the series using
|
|
Packit |
67cb25 |
a Levin :math:`u`-transform. Additional working space must be provided in
|
|
Packit |
67cb25 |
:data:`w`. The extrapolated sum is stored in :data:`sum_accel`. The actual
|
|
Packit |
67cb25 |
term-by-term sum is returned in :code:`w->sum_plain`. The algorithm
|
|
Packit |
67cb25 |
terminates when the difference between two successive extrapolations
|
|
Packit |
67cb25 |
reaches a minimum or is sufficiently small. The difference between these
|
|
Packit |
67cb25 |
two values is used as estimate of the error and is stored in
|
|
Packit |
67cb25 |
:data:`abserr_trunc`. To improve the reliability of the algorithm the
|
|
Packit |
67cb25 |
extrapolated values are replaced by moving averages when calculating the
|
|
Packit |
67cb25 |
truncation error, smoothing out any fluctuations.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
Examples
|
|
Packit |
67cb25 |
========
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
The following code calculates an estimate of :math:`\zeta(2) = \pi^2 / 6`
|
|
Packit |
67cb25 |
using the series,
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. math:: \zeta(2) = 1 + 1/2^2 + 1/3^2 + 1/4^2 + \dots
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
After :data:`N` terms the error in the sum is :math:`O(1/N)`, making direct
|
|
Packit |
67cb25 |
summation of the series converge slowly.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. include:: examples/sum.c
|
|
Packit |
67cb25 |
:code:
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
The output below shows that the Levin :math:`u`-transform is able to obtain an
|
|
Packit |
67cb25 |
estimate of the sum to 1 part in
|
|
Packit |
67cb25 |
:math:`10^{10}`
|
|
Packit |
67cb25 |
using the first eleven terms of the series. The
|
|
Packit |
67cb25 |
error estimate returned by the function is also accurate, giving
|
|
Packit |
67cb25 |
the correct number of significant digits.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. include:: examples/sum.txt
|
|
Packit |
67cb25 |
:code:
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
Note that a direct summation of this series would require
|
|
Packit |
67cb25 |
:math:`10^{10}`
|
|
Packit |
67cb25 |
terms to achieve the same precision as the accelerated
|
|
Packit |
67cb25 |
sum does in 13 terms.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
References and Further Reading
|
|
Packit |
67cb25 |
==============================
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
The algorithms used by these functions are described in the following papers,
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
* T. Fessler, W.F. Ford, D.A. Smith,
|
|
Packit |
67cb25 |
HURRY: An acceleration algorithm for scalar sequences and series
|
|
Packit |
67cb25 |
*ACM Transactions on Mathematical Software*, 9(3):346--354, 1983.
|
|
Packit |
67cb25 |
and Algorithm 602 9(3):355--357, 1983.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
The theory of the :math:`u`-transform was presented by Levin,
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
* D. Levin,
|
|
Packit |
67cb25 |
Development of Non-Linear Transformations for Improving Convergence of
|
|
Packit |
67cb25 |
Sequences, *Intern.: J.: Computer Math.* B3:371--388, 1973.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
A review paper on the Levin Transform is available online,
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
* Herbert H. H. Homeier, Scalar Levin-Type Sequence Transformations,
|
|
Packit |
67cb25 |
http://arxiv.org/abs/math/0005209
|