Package levin provides functions that implement the Levin u-transformation (Levin 1973). This transformation may accelerate convergence of the sum of infinite sequences using a small number of terms.
Methods to accelerate the convergence of infinite series do not usually return the exact sum. There are two major sources of error: truncation error, due to the finite number of terms used; and, rounding error, due to the finite precision arithmetic used in the calculations. Truncation error dominates for a small number of terms (perhaps 5 to 10 terms). Rounding error dominates as the number of terms increases (perhaps above 10 to 20 terms for IEEE-754 double precision arithmetic).
This package has two functions that address these sources of error.
Function levin_u_sum (by default) uses rational arithmetic with a
fixed number
of terms. This eliminates rounding error if the terms are known exactly.
Function bflevin_u_sum uses bigfloat precision.
It uses error estimates to increase the number of terms and the bigfloat
arithmetic precision to meet an accuracy target.
See Examples for levin for examples.
See References for levin for further information.
Estimate sum(a(n), n, n_0, inf) using at most nterms
terms using the Levin u-transform see Levin 1973.
The following values are recognized for the optional argument mode.
If mode is not supplied, it is assumed to be levin_algebraic.
levin_algebraicThe calculation is performed in exact arithmetic. levin_u_sum returns the result.
levin_numericThe calculation is performed in bigfloat arithmetic.
The return value is a list [result, variance]
where result is the result of the bigfloat calculation,
and variance is in units of 10^(-2*fpprec).
load("levin") loads this function.
See Examples for levin for examples.
Estimate sum(a(n), n, n_0, inf) using the
Levin u-transform in bigfloat arithmetic.
bflevin_u_sum attempts to return the sum of the infinite series
with a precision given by the global variable fpprec using bigfloat arithmetic.
See levin_options for options to control this function.
bflevin_u_sum uses an adaptive algorithm to increase both the number of
terms and the bigfloat precision used for internal calculations
until the estimated error is acceptable.
load("levin") loads this function.
See Examples for levin for examples.
Function bflevin_u_sum attempts to return the sum of an infinite series
with a precision given by the global variable fpprec using bigfloat arithmetic.
bflevin_u_sum uses an adaptive algorithm to increase both the number of
terms used and the bigfloat precision used for internal calculations
until the estimated error is acceptable.
The undeclared array levin_options contains options for controlling bflevin_u_sum.
Note that the subscript values for levin_options are strings.
levin_options["debug"]When true, bflevin_u_sum generates additional output. Default: false
levin_options["min_terms"]Minimum number of terms used by bflevin_u_sum. Default: 5
levin_options["max_terms"]Maximum number of terms used by bflevin_u_sum. Default: 640 (equal to 5*2^7)
levin_options["min_precision"]Initial bigfloat precision for bflevin_u_sum. Default: 16
levin_options["max_precision"]Maximum bigfloat precision for bflevin_u_sum. Default: 1000
The Basel problem asks for the precise summation of the reciprocals of the squares of the natural numbers. Leonhard Euler found, in 1734,
\[{\sum_{n=1}^{\infty}{{1}\over{n^2}}}
= {{\pi^2}\over{6}}
= {1.6449340668482262...}
\]
The sum is evaluated using both: levin_u_sum using modes levin_algebraic and levin_numeric; and
bflevin_u_sum with two values of bigfloat precision fpprec.
(%i1) (load("levin"), exact: %pi^2/6, done);
(%o1) done
(%i2) s: levin_u_sum(1/n^2, n, 1, 10);
3899836039
(%o2) ----------
2370816000
(%i3) float(s); (%o3) 1.644934081345832
(%i4) float(s - exact); (%o4) 1.4497605782537448e-8
(%i5) s: levin_u_sum(1/n^2, n, 1, 10, 'levin_numeric); (%o5) [1.644934081345756b0, 1.882785043290232b-12]
(%i6) bfloat(s[1] - exact); (%o6) 1.449752928817105b-8
(%i7) s: bflevin_u_sum(1/n^2, n, 1); (%o7) 1.644934066848226b0
(%i8) s - bfloat(exact); (%o8) 2.775557561562891b-17
(%i9) /* Now increase fpprec and try the same example again. */ fpprec: 32; (%o9) 32
(%i10) s: bflevin_u_sum(1/n^2, n, 1); (%o10) 1.644934066848226436472415166646b0
(%i11) s - bfloat(exact); (%o11) - 3.0814879110195773648895647081359b-33
Using 10 terms in the series we were able to extrapolate to the exact sum with an error of approximately \(10^{-8}\). This would require around \(10^{8}\) terms by direct summation. In this case the two modes of levin_u_sum returned results of similar accuracy.
The effect of the number of terms nterms on the accuracy of levin_u_sum is shown in the following example. The sum of the series and the approximation error is evaluated for increasing values of argument nterms for both values of the optional argument mode: levin_algebraic and levin_numeric. The numeric calcuations are performed with bigfloat precision fpprec of 16. Errors are calculated with fpprec equal to 64.
The results are reported in three columns: nterms, the number of terms used; errora, the difference between the algebraic result and the exact sum; and errorn, the difference between numeric result and the exact sum. For small values of nterms the two modes return similar results. As nterms increases above 12, the error for mode levin_algebraic continues to decrease. However, for mode levin_numeric the error increases with nterms to due to roundoff. This behavior is commonly encountered, and can be addressed by increasing fpprec.
(%i1) (load("levin"), fpprec: 32, fpprintprec: 3, exact: %pi^2/6, done);
(%o1) done
(%i2) for nterms: 6 step 2 thru 24 do (
sa: levin_u_sum(1/n^2, n, 1, nterms),
sn: block([fpprec: 16], levin_u_sum(1/n^2, n, 1, nterms, 'levin_numeric)),
errora: abs(bfloat(sa - exact)),
errorn: abs(bfloat(sn[1] - exact)),
print(nterms, " ", errora, " ", errorn));
6 2.58b-4 2.58b-4
8 3.51b-6 3.51b-6
10 1.45b-8 1.45b-8
12 2.44b-10 2.44b-10
14 5.16b-12 3.78b-11
16 3.54b-14 4.5b-10
18 3.75b-16 3.21b-11
20 1.36b-17 8.66b-9
22 1.73b-19 4.48b-7
24 4.45b-22 4.25b-6
(%o2) done
The Levin u-transformation can find the anti-limit of certain divergent series. In this example we calculate the anti-limit of the divergent Taylor series for log(1+x) when x=2.
\[{A_n} = {-\sum_{m=1}^{n}{{\left( -2 \right)}^{m}\over{m}}}
\]
with anti-limit equal to log(3) = 1.0986122886681098...
(%i1) (load("levin"), done);
(%o1) done
(%i2) levin_u_sum(-(-2)^n/n, n, 1, 10);
3641179708
(%o2) ----------
3314344635
(%i3) s: float(%); (%o3) 1.0986122775370342
(%i4) exact: log(3.0); (%o4) 1.0986122886681098
(%i5) s - exact; (%o5) - 1.1131075616788166e-8
The t, u and v Levin transformations are introduced in:
Smith and Ford compared the performance of series acceleration algorithms in a series of papers. In conjunction with Fessler they published algorithms and code that implement the u-transformation with estimates of truncation and round-off errors.
Sidi (2003) discusses u transformations and more recent extensions.
Bender and Orszag (1978), Chapter 8 Summation of Series, provides general background on motivation and methods, but does not cover Levin transformations.