Numerical evaluation of a sum

Question

I am trying to compute numerically NSum[(-1)^n/n^3, {n, 1, Infinity}]. Of course, using first Sum would work here, but often it's not an option (no closed-form). I would like to compute the sum with an arbitrary number of digits.

So far, I've found that

NSum[(-1)^n/n^3, {n, 1, Infinity},
 WorkingPrecision -> 24, 
 AccuracyGoal -> 24, PrecisionGoal -> Infinity]

will give only 16 digits (-0.9015426773696957), whatever the value of WorkingPrecision and AccuracyGoal.

However, if I add the option Method -> "AlternatingSigns", I get as many digits as I want, and it's really controlled by WorkingPrecision.

Among the other options, Method -> "EulerMaclaurin" gives a wrong answer. Not all digits are correct, and there is an imaginary part: -0.901542677369695714049804 + 15.503138340149910087738157 I. And Method -> "WynnEpsilon" also fails, giving -0.9015426773696957140, even if I ask for more digits.

So, is there something I'm doing wrong? Is there a reason it does not work if I don't give a specific method? Is there a more straightforward way to compute a sum (or any numerical result, for that matter)?

And, since I'm quite interested by numerical computations with Mathematica and I'm a beginner (with Mathematica, not with numerical analysis :-)), is there a good guide out there ?

As an aside, I know Maxima better, and I'm used to the option fpprec, that controls "bigfloat" evaluation in every computations. Is there something similar in Mathematica? I think it has a better control of accuracy, so maybe the question does not make sense (maybe there is no "global" setting).

---

To answer ciao's comment, since there is not much space in comments:

Actually, if I use larger values of the option NSumTerms, I get more and more digits, but they don't agree with what I get with Method -> "AlternatingSigns":

NSum[(-1)^n/n^3, {n, 1, Infinity}, WorkingPrecision -> 50, 
 AccuracyGoal -> 50, PrecisionGoal -> Infinity, NSumTerms -> 5000]

-0.90154267736969571404980362113358749307373971926

NSum[(-1)^n/n^3, {n, 1, Infinity}, WorkingPrecision -> 50, 
 AccuracyGoal -> 50, PrecisionGoal -> Infinity, 
 Method -> "AlternatingSigns"]

-0.90154267736969378050491347983360150530995076960417

After 14 digits, they are different numbers.

And here it's possible to check:

Sum[(-1)^n/n^3, {n, 1, Infinity}] // N[#, 50] &

-0.90154267736969571404980362113358749307373971925537

So, the answer with Method -> "AlternatingSigns" is also wrong!

And it's not better if I combine Method and NSumTerms options:

NSum[(-1)^n/n^3, {n, 1, Infinity}, WorkingPrecision -> 50, 
 AccuracyGoal -> 50, PrecisionGoal -> Infinity, 
 Method -> "AlternatingSigns", NSumTerms -> 5000]

-0.90154267736969571404546474460914681065800480608588

So the results differ after 20 digits.

Now I'm even more puzzled :-) I can see that your comment is right, since I get the right answer, but I absolutely don't understand why, and why the other methods fail without notification.

Also, it's not a very satisfying approach, since I can't set the number of accurate digits in the output (it depends on NSumTerms in a non obvious way).

Answer 1

As a general rule, one should usually set WorkingPrecision to be higher than either of AccuracyGoal and PrecisionGoal (if not set to ∞) to mitigate any roundoff error incurred during the evaluation.

For the specific methods implemented within NSum[]: it is documented, but not apparently well-known, that you can change the way the integral term in the Euler-Maclaurin method is evaluated numerically. For instance,

NSum[(-1)^n/n^3, {n, 1, ∞},
     WorkingPrecision -> 30, AccuracyGoal -> 24, PrecisionGoal -> ∞, 
     Method -> {"EulerMaclaurin",
                Method -> {NIntegrate, Method -> "GlobalAdaptive"}}]
   -0.901542677369695714049804 + 0.*10^-25 I

yields a result quite close to $\frac34\zeta(3)$.

For the Wynn $\varepsilon$ method, what you can tweak is the degree of the transformation used (the default Degree -> 1 corresponds to the Aitken $\Delta^2$ process) and the number of terms of the sum to be used in the extrapolation, "ExtraTerms". (Recall that in evaluating a sum by extrapolation, one usually splits the sum in two parts, where one part is evaluated through direct summation, and the "tail" is evaluated by using extrapolation.)

Thus,

NSum[(-1)^n/n^3, {n, 1, ∞},
     WorkingPrecision -> 30, AccuracyGoal -> 24, PrecisionGoal -> ∞,
     NSumTerms -> 20, Method -> {"WynnEpsilon", Degree -> 2, "ExtraTerms" -> 20}]
   -0.901542677369695714049803590206

Of course, NSum[] is not always able to return good answers, so anybody who is sufficiently paranoid about his results should always look into a number of other summation methods. But that is for another time for me to expound on...

Answer 2

Comments

The author of the question says this in a comment:

Is there a way to simply tell Mathematica "I want 50 digits, choose the options as you wish"? Maybe I'm asking for too much, but I find rather strange that many options give wrong results, and one option give the right answer, but with the correct parameter.

Here are some comments of mine in that regard.

1. The Numerical Summation Extrapolation (NSE) methods seem to have narrow scopes of applicability within the "space" of number sequences. If the model of a NSE method fits the sequence that we work with, then the NSE method is usually very effective. That same method can be failing really badly with sequences that do not fit its model.

2. From that perspective NSum's set of methods is deficient. E.g. NSum does not have Richardson extrapolation.

3. In general Mathematica tries to provide algorithms that select the methods and parameters to provide the answers the users want. So with "I want 50 digits, choose the options as you wish" you are not asking too much from Mathematica (generally speaking) but you are asking too much from NSum.

Code

Here is an implementation that solves the particular problem posted in the question. It uses the Shanks transformation, and a search procedure that is using the difference of consecutive results for a stopping criteria. (As in Cauchy's convergence test.) I tried the Richardson extrapolation too but it did not give good results. (Using the Shanks-vs-Richardson comparison code in this answer of mine.)

First we define the Shanks transformation function:

ClearAll[Shanks]
Shanks[A_, n_] := (A[n + 2] A[n] - A[n + 1]^2)/(A[n + 2] + A[n] - 2 A[n + 1]);
Shanks[A_, 1, n_] := Shanks[A, 1, n] = Shanks[A, n];
Shanks[A_, k_, n_] := 
  Shanks[A, k, n] = ((Shanks[A, -1 + k, n] Shanks[A, -1 + k, 2 + n] - 
       Shanks[A, -1 + k, 1 + n]^2)/(Shanks[A, -1 + k, n] + 
       Shanks[A, -1 + k, 2 + n] - 2 Shanks[A, -1 + k, 1 + n]));

Here is the search function:

Clear[SearchSumValue]
SearchSumValue[partialSumFunc_, accGoal_, step_: 10, startStep_: 40, methodFunc_: Shanks] :=
  Block[{res, pf = 2},
   res =
    NestWhile[
     {#[[1]] + step, #[[3]], 
       methodFunc[partialSumFunc, #[[1]] + step]} &,
     {1, methodFunc[partialSumFunc, startStep], 
      methodFunc[partialSumFunc, startStep + step]},
     Abs[N[#[[2]], pf*accGoal] - N[#[[3]], pf*accGoal]] > 
       10^-accGoal &];
   (* n-steps, estimate, error *)
   Append[res[[{1, 3}]], 
    Abs[N[res[[2]], pf*accGoal] - N[res[[3]], pf*accGoal]]]
   ];

Here is the application of the functions above:

sf[n_] := Sum[(-1)^i/i^3, {i, 1, n}]

(* partialSumFunc, accGoal, step, startStep, methodFunc *)
res = SearchSumValue[sf, 20, 30, 40, Shanks[#1, 6, #2] &];
N[res[[2]] - exact, 40]  
(* 8.757707518816394163494370352471694535835*10^-28 *)

(* partialSumFunc, accGoal, step, startStep, methodFunc *)
res = SearchSumValue[sf, 50, 40, 40, Shanks[#1, 6, #2] &];
N[res[[2]] - exact, 50]
(* -3.1179069154714173193769908669678269443598103298019*10^-50 *)

Answer 3

x = N[Sum[(-1)^n/n^3, {n, 1, Infinity}], 47]


y = NSum[(-1)^n/n^3, {n, 1, Infinity}, WorkingPrecision -> 47, 
  AccuracyGoal -> 47, PrecisionGoal -> 47, NSumTerms -> 10000]

z = NSum[(-1)^n/n^3, {n, 1, Infinity}, Method -> "AlternatingSigns", 
  WorkingPrecision -> 47, AccuracyGoal -> 47, PrecisionGoal -> 47]

x == y == z

True (-0.90154267736969571404980362113358749307373971926)

Using the simplest form:

e = NSum[(-1)^n/n^3, {n, 1, Infinity}]

we get with

e // FullForm

-0.9015426773696957`

which is clearly less precise. But!

x === y === z === e

True

Strange ...

Answer 4

This doesn't directly answer the original question, but why not use the fact that Mathematica can evaluate the exact sum symbolically:

    s = Sum[(-1)^n/n^3, {n, 1, Infinity}] 
(-3*Zeta[3])/4 *)

Then you may use N, with a 2nd argument to specify the desired precision, e.g.:

    N[s, 50]
(* -0.90154267736969571404980362113358749307373971925537 *)

(Of course such an approach will only work with those special series whose sums Mathematica already knows about somehow.)

Answer 5

Edit: [I see that someone else gave you the same information I give in the following, about NSumTerms. I hope the following will at least reinforce the solution in your work.]

I found NSumTerms very helpful when Mathmatica becomes intransigent in my research, which is often! NSumTerms works well on your problem too. AlternatingSigns usually is necessary on stubborn alternating series that have no or very complicated closed forms like,

Sum[(1)^n (n^(1/n) - 1), {n, 1, Infinity}] 

ie.

NSum[(-1)^x (x^(1/x) - 1), {x, 1, Infinity}, Method -> "AlternatingSigns", WorkingPrecision -> 80]

or

Sum[(-1)^x ((Log[x]/x)^n/n!), {x, 1, Infinity}, {n, 1, Infinity}] 

ie.

 NSum[(-1)^x (Log[x]/x)^n/n!, {x, 1, Infinity}, {n, 1, Infinity}, 
  Method -> "AlternatingSigns", WorkingPrecision -> 80, 
  PrecisionGoal -> 80, AccuracyGoal -> 80, NSumTerms -> 80]

. Both of the just mentioned codes give rise to approximations of the MRB constant, the most fruitful result of my research.

In your case, you need a whoping 7'th power of 10 terms, or there about, to get 80 accurate digits, like in the following:

NSum[(-1)^n/n^3, {n, 1, Infinity}, NSumTerms -> 10^7, 
 WorkingPrecision -> 80, AccuracyGoal -> 80, 
 PrecisionGoal -> Infinity]

gives

-0.9015426773696957140498036211335874930737397192553741613442036665063\ 7865433973482.

I am sure it is correct because if you compute the digits from the closed form, (as given above), like in the following...

N[-3/4 Zeta[3], 80]

you also get

-0.9015426773696957140498036211335874930737397192553741613442036665063\ 7865433973482

Alternating signs actually gives only 31 digits of accuracy in the following attempt.

 NSum[(-1)^n/n^3, {n, 1, Infinity}, NSumTerms -> 10^7, 
 Method -> "AlternatingSigns", WorkingPrecision -> 80, 
 AccuracyGoal -> 80, PrecisionGoal -> Infinity]

Answer 6

I know you want to use NSum but there are other options for accelerating the convergence of series.

Start by noting that your series

$\left\{-1,\frac{1}{8},-\frac{1}{27},\frac{1}{64},-\frac{1}{125},\frac{1}{216},-\frac{1}{343},...\right\}$

is the negative of this one

$\left\{1,-\frac{1}{8},\frac{1}{27},-\frac{1}{64},\frac{1}{125},-\frac{1}{216},\frac{1}{343},...\right\}$

which has the general term

$\frac{(-1)^k}{(k+1)^3}$

for k = 0, 1, 2,...

There is a little routine I came across which I thought was from the Borweins (see J.M.'s comment below), that I call the converger.

acc[n_]:=Module[{d,b,c,s},
d=(3+Sqrt[8])^n;
d=(d+1/d)/2;
b=-1;
c=-d;
s=0;
Table[c=b-c;s=s+c*a[k];b=(k+n)(k-n) b/((k+1/2)(k+1)),{k,0,n-1}];
s/d]

You prime it by entering the general term as a positive sequence:

a[k_]:=1/(k+1)^3

Now just run

-N[acc[300],200]

and you have about 200 places. The parameter n controls the number of iterations and hopefully the convergence. Quite fast and might get you out of jam when nothing else is working.

Answer 7

If you look at the FullForm or InputForm of the NSum result without using "AlternatingSigns" you will see that more digits are available; however, display is limited by the determined precision which is less than 17 digits. The determined precision which is implied by the shortened result also can be seen using Precision.

approx1 = NSum[(-1)^n/n^3, {n, 1, Infinity}, WorkingPrecision -> 24, AccuracyGoal -> 24, PrecisionGoal -> Infinity]

-0.9015426773696957

#[approx1] & /@ {Precision, FullForm, InputForm}

{16.774,-0.9015426773696956967432522482628641510116.773986512910163,-0.9015426773696956967432522482628641510116.773986512910163}

If the precision of the result is less than the requested precision then you should start investigating other methods.

approx2 = NSum[(-1)^n/n^3, {n, 1, Infinity}, WorkingPrecision -> 24, AccuracyGoal -> 24, PrecisionGoal -> Infinity, Method -> "AlternatingSigns"]

-0.901542677369693780504913

#[approx2] & /@ {Precision, FullForm, InputForm}

{24.1487,-0.90154267736969378050491347983360150530995076960472861836724.148671096657576,-0.90154267736969378050491347983360150530995076960472861836724.148671096657576}