Function Repository Resource:

LevinSum

Source Notebook

Evaluate an infinite sum using the Levin transformation

Contributed by: Jan Mangaldan

ResourceFunction["LevinSum"][f,{i,imin,}]

numerically evaluates the sum using the Levin transformation.

Details and Options

Levin's transformation is based on the evaluation of the expression , where Sn is the nth partial sum of the series, gn is an auxiliary sequence representing an error estimate of the infinite sum, and is DifferenceDelta[f,{i,n}].
The following options can be given:
"ExtraTerms"15number of terms to use in the Levin transform
"Terms"15number of terms to sum directly
"Type"Automaticthe type of Levin transformation to use
WorkingPrecisionMachinePrecisionthe precision used in internal computations
"Type" can be set to any one of the following types, corresponding to different Levin transforms:
"T"t-transformation, gn=Sn-Sn-1
"U"u-transformation, gn=(n+1)(Sn-Sn-1)
"V"v-transformation, gn=-(Sn+1-Sn)(Sn-Sn-1)/(Sn+1-2Sn+Sn-1)
"D"d-transformation, gn=Sn+1-Sn
The t and d transformations work best with alternating series, while the u and v transformations are more flexible.
You should realize that with sufficiently pathological summands, ResourceFunction["LevinSum"] can give wrong answers. In most cases, you can test the answer by looking at its sensitivity to changes in the setting of options for ResourceFunction["LevinSum"].
ResourceFunction["LevinSum"] has attribute HoldAll, and effectively uses Block to localize variables.

Examples

Basic Examples (2) 

Evaluate the alternating harmonic series:

In[1]:=
ResourceFunction["LevinSum"][(-1)^(j - 1)/j, {j, 1, \[Infinity]}]
Out[1]=

Compare with the closed form:

In[2]:=
% - Log[2]
Out[2]=

Options (8) 

ExtraTerms (2) 

Use 25 terms for the Levin transformation:

In[3]:=
ResourceFunction["LevinSum"][(-1)^k/(2 k + 1), {k, 0, \[Infinity]}, "ExtraTerms" -> 25, WorkingPrecision -> 25]
Out[3]=

Compare with the exact result:

In[4]:=
% - \[Pi]/4
Out[4]=

Set "Terms" to 0 so that all terms are used in extrapolation:

In[5]:=
ResourceFunction["LevinSum"][(-1)^k/(2 k + 1), {k, 0, \[Infinity]}, "ExtraTerms" -> 25, "Terms" -> 0, WorkingPrecision -> 35]
Out[5]=

Compare with the exact result:

In[6]:=
% - \[Pi]/4
Out[6]=

Terms (2) 

Directly sum the first 25 terms before applying the Levin transformation:

In[7]:=
ResourceFunction["LevinSum"][(-1)^k/(2 k + 1), {k, 0, \[Infinity]}, "Terms" -> 25, WorkingPrecision -> 25]
Out[7]=

Compare with the exact result:

In[8]:=
% - \[Pi]/4
Out[8]=

Type (2) 

Show the results of the different Levin transformations on an alternating series:

In[9]:=
TableForm[
 Table[With[{r = ResourceFunction["LevinSum"][(-1)^k/(
      2 k + 1), {k, 0, \[Infinity]}, "Type" -> t, WorkingPrecision -> 25]}, {t, r, r - \[Pi]/4}], {t, {"T", "U", "V", "D"}}], TableHeadings -> {None, {"Type", "Result", "Error"}}]
Out[9]=

Show the results of the different Levin transformations on a non-alternating series:

In[10]:=
TableForm[
 Table[With[{r = ResourceFunction["LevinSum"][1/k^2, {k, 1, \[Infinity]}, "Type" -> t, WorkingPrecision -> 25]}, {t, r, r - \[Pi]^2/6}], {t, {"T", "U", "V", "D"}}], TableHeadings -> {None, {"Type", "Result", "Error"}}]
Out[10]=

WorkingPrecision (2) 

Use a higher setting of WorkingPrecision:

In[11]:=
ResourceFunction["LevinSum"][(-1)^k/(2 k + 1), {k, 0, \[Infinity]}, WorkingPrecision -> 25]
Out[11]=

Compare with the exact result:

In[12]:=
% - \[Pi]/4
Out[12]=

Applications (2) 

Use the Levin d-transform to evaluate the Dirichlet eta function:

In[13]:=
eta[z_?NumericQ] := ResourceFunction["LevinSum"][(-1)^(j - 1)/j^z, {j, 1, \[Infinity]}, "Type" -> "D"]

Compare with the built-in DirichletEta:

In[14]:=
Plot[{DirichletEta[z], eta[z]}, {z, -3, 3}, PlotStyle -> {Thick, Dashed}]
Out[14]=

Plot the relative error:

In[15]:=
Plot[Abs[DirichletEta[z] - eta[z]]/Abs[DirichletEta[z]], {z, -3, 3}, PlotRange -> All]
Out[15]=

Use the Levin v-transform with NIntegrate to numerically evaluate an oscillatory integral:

In[16]:=
term[j_Integer] := NIntegrate[
  BesselJ[0, x], {x, BesselJZero[0, j], BesselJZero[0, j + 1]}, WorkingPrecision -> 25]
NIntegrate[BesselJ[0, x], {x, 0, BesselJZero[0, 1]}, WorkingPrecision -> 25] + ResourceFunction["LevinSum"][term[j], {j, 1, \[Infinity]}, "Type" -> "V", WorkingPrecision -> 25]
Out[16]=

Compare with the exact result:

In[17]:=
% - 1
Out[17]=

Properties and Relations (2) 

Directly summing the first few terms of a series usually does not give sufficient accuracy:

In[18]:=
Sum[N[(-1)^k/(2 k + 1), 35], {k, 0, 19}]
Out[18]=
In[19]:=
% - \[Pi]/4
Out[19]=

Using the Levin transform on a series often gives better results:

In[20]:=
ResourceFunction["LevinSum"][(-1)^k/(2 k + 1), {k, 0, \[Infinity]}, "ExtraTerms" -> 20, "Terms" -> 0, WorkingPrecision -> 35]
Out[20]=
In[21]:=
% - \[Pi]/4
Out[21]=

Possible Issues (2) 

LevinSum may give finite results for formally divergent series:

In[22]:=
ResourceFunction["LevinSum"][(-1)^j j!, {j, 0, \[Infinity]}, "Terms" -> 0, WorkingPrecision -> 25]
Out[22]=

Compare with the exact result:

In[23]:=
% - E ExpIntegralE[1, 1]
Out[23]=

Neat Examples (2) 

Numerically evaluate a formally divergent oscillatory integral:

In[24]:=
est = ResourceFunction["LevinSum"][
   NIntegrate[Sin[x] Log[x], {x, j \[Pi], (j + 1) \[Pi]}, WorkingPrecision -> 25], {j, 0, \[Infinity]}, "Type" -> "D", WorkingPrecision -> 25] // Quiet
Out[24]=

Compare with the exact answer:

In[25]:=
ex = Limit[\!\(
\*SubsuperscriptBox[\(\[Integral]\), \(0\), \(\[Infinity]\)]\(
\*SuperscriptBox[\(E\), \(\(-\[CurlyEpsilon]\)\ x\)] Sin[x] Log[
     x] \[DifferentialD]x\)\), \[CurlyEpsilon] -> 0, Direction -> "FromAbove"]
Out[25]=
In[26]:=
est - ex
Out[26]=

Version History

  • 1.0.0 – 05 April 2021

Source Metadata

Related Resources

License Information