8
$\begingroup$

I'm doing a sanity check of the following equation: $$\sum_{j=2}^\infty \frac{(-x)^j}{j!}\zeta(j) \approx x(\log x + 2 \gamma -1)$$

Naive comparison of the two shows a bad match but I suspect one of the graphs is incorrect.

  1. Why isn't there a warning?
  2. How do I compute this sum correctly?
katsurda[x_] := NSum[(-x)^j/j!  Zeta[j], {j, 2, Infinity}];
katsurdaApprox[x_] := x (Log[x] + 2 EulerGamma - 1) - Zeta[0];
plot1 = DiscretePlot[katsurda[x], {x, 0, 40, 2}];
plot2 = Plot[katsurdaApprox[x], {x, 0, 40}];
Show[plot1, plot2]

enter image description here

  1. meta How do I avoid being mislead by incorrect numeric results? Would using NIntegrate instead of NSum give better guarantees? My usual approach of a avoiding machine precision, checking Precision of the answer and minding warnings fails in the example below
katsurda[x_] := 
  NSum[(-x)^j/j! Zeta[j], {j, 2, Infinity}, WorkingPrecision -> 32, 
   NSumTerms -> 2.5 x];
katsurdaApprox[x_] := x (Log[x] + 2 EulerGamma - 1) - Zeta[0];
Print["Precision: ", Precision@katsurda[100]] (* 13.9729 *)
Print["Discrepancy: ", katsurda[100] - katsurdaApprox[100]] (* 94.65088290385, but should be <1 *) 

Background: the expression comes from "Power series with the Riemann zeta-function in the coefficients" by Katsurada M (paper)

$\endgroup$
4
  • 1
    $\begingroup$ BTW, I don't know a general answer for question 3 for NSum. I have a good idea how many of the methods of NIntegrate and NDSolve work. One test is to raise WorkingPrecision and PrecisionGoal, and perhaps AccuracyGoal, and see if the result is stable. That hasn't always been reliable for NSum for me. The Precision of an answer in Mma usually indicates only the effect of round-off. It is not good indicator of convergence/nonconvergence. For the sum at hand, since it is alternating and eventually decreasing in magnitude, Method -> "AlternatingSigns", should be reliable. $\endgroup$
    – Michael E2
    Jul 8, 2021 at 16:25
  • $\begingroup$ Perhaps the issue is the internal heuristics governed by NSumTerms. I tried to use NIntegrate by splitting series into positive/negative parts and treating each as Riemann sum, but this approximation turned out to be bad for large x wolframcloud.com/obj/yaroslavvb/newton/forum-katsurda-split.nb $\endgroup$ Jul 8, 2021 at 16:33
  • $\begingroup$ I did something similar, but using successive pairs to get a non-alternating series. NIntegrate couldn't handle the it. $\endgroup$
    – Michael E2
    Jul 8, 2021 at 16:42
  • $\begingroup$ @Daniel Lichtblau you guys should use this to debug NSumTerms->Automatic :) $\endgroup$ Jul 8, 2021 at 16:45

4 Answers 4

9
$\begingroup$

The sum is alternating, so you might need extra precision and NSumTerms:

katsurda[x_] := 
  NSum[(-x)^j/j! Zeta[j], {j, 2, Infinity}, WorkingPrecision -> 16, 
   NSumTerms -> Max[15, 2 x]];
katsurdaApprox[x_] := x (Log[x] + 2 EulerGamma - 1) - Zeta[0];
plot1 = DiscretePlot[katsurda[x], {x, 0, 40, 2}];
plot2 = Plot[katsurdaApprox[x], {x, 0, 40}];
Show[plot1, plot2, PlotRange -> All]

enter image description here

Note: There is a warning, but it's suppressed by the plotter:

Block[{x = 30},
 NSum[(-x)^j/j! Zeta[j], {j, 2, Infinity}]
 ]

NumericalMath`NSequenceLimit::seqlim: The general form of the sequence could not be determined, and the result may be incorrect.

(*  126.442  *)

Update:

Here's a way to estimate the needed precision by estimating the largest term in the series. The method "AlternatingSigns" is (or should be) a reliable method for the sum, provided the working precision is high enough.

katsurda[x_] := NSum[(-x)^j/j! Zeta[j], {j, 2, Infinity},
   WorkingPrecision -> 
    16 + 
     FindMaxValue[{(j0*Log[x] - LogGamma[1 + j0])/Log[10], 
       j0 > 1}, {j0, x}], NSumTerms -> Max[15, 4 + 2 x],
   Method -> "AlternatingSigns"];
katsurdaApprox[x_] := x (Log[x] + 2 EulerGamma - 1) - Zeta[0];
plot1 = DiscretePlot[katsurda[x], {x, Round[1.5^Range[10, 17]]}];
plot2 = Plot[katsurdaApprox[x], {x, 50, 1000}];
Show[plot1, plot2, PlotRange -> All]

enter image description here

Update 2:

After some playing, we can see that the maximum for large x is around j == x, and Method -> "AlternatingSigns" adjusts the number of terms needed automatically. So WorkingPrecision -> 16 + x ensures a sufficient precision for x >= 0. Thus here's a simplified code:

katsurda[x_] := 
 NSum[(-x)^j/j! Zeta[j], {j, 2, Infinity}, 
  Method -> "AlternatingSigns", WorkingPrecision -> 16 + x]
$\endgroup$
4
  • $\begingroup$ Interesting! BTW, the number $16$ seems special here, how did you come up with it? If I use WorkingPrecision -> 32 then I get a big discrepancy and no warning in katsurda[100] - katsurdaApprox[100] $\endgroup$ Jul 7, 2021 at 23:25
  • 1
    $\begingroup$ @YaroslavBulatov I'm just figuring out that 2 x for NSumTerms is not large enough for large x. I'm not sure how to determine the "right" number. At least yet :) -- Sorry, 2.5 x doesn't really work for x = 100. $\endgroup$
    – Michael E2
    Jul 7, 2021 at 23:27
  • $\begingroup$ Using NSumTerms -> 2.5 x terms and WorkingPrecision->32 tells me that this approximation is off by 94.65 for $x=100$. There's no warning and 94.65 answer has 13.97 digits of precision. But I think this answer is wrong....lack of warning bothers me a bit... $\endgroup$ Jul 7, 2021 at 23:35
  • $\begingroup$ @YaroslavBulatov Yes, I had to use WorkingPrecision -> 64 Or 128? Clearly both NSumTerms and WorkingPrecision have to go up with x. Also Method -> "AlternatingSigns", which it does not default to, seems to work better (and maybe faster). $\endgroup$
    – Michael E2
    Jul 7, 2021 at 23:41
2
$\begingroup$

Using the $d$-type Weniger transformation, as implemented in ResourceFunction["WenigerSum"], we get the following:

Plot[{ResourceFunction["WenigerSum"][(-x)^j Zeta[j]/j!, {j, 2, ∞},
                                     "ExtraTerms" -> 25, "Type" -> "D", WorkingPrecision -> 20], 
      x (Log[x] + 2 EulerGamma - 1) - Zeta[0]} // Evaluate, {x, 0, 45}, WorkingPrecision -> 25]

Weniger sum vs. asymptotic approximation

Plot[ResourceFunction["WenigerSum"][(-x)^j Zeta[j]/j!, {j, 2, ∞},
                                    "ExtraTerms" -> 25, "Type" -> "D", WorkingPrecision -> 20] -
     (x (Log[x] + 2 EulerGamma - 1) - Zeta[0]) // Evaluate, {x, 0, 45},
     PlotRange -> All, WorkingPrecision -> 25]

residual

One would need to play around with the "ExtraTerms" and "Terms" settings for this function (e.g. try the combination "Terms" -> 0, "ExtraTerms" -> 35), and one has to use high precision as well to stave off the looming threat of catastrophic cancellation, especially for large $x$. I have chosen the $d$-type transformation since the original series is alternating (as recommended by Weniger), but the $t$-type ("Type" -> "T") can be used on such series as well. I invite you to do your own experiments.

$\endgroup$
1
$\begingroup$

Turn the sum around to make it non-alternating:

$$ \sum_{j=2}^{\infty}\frac{(-x)^j}{j!}\zeta(j) = \sum_{j=2}^{\infty}\frac{(-x)^j}{j!}\left(\sum_{n=1}^{\infty}\frac{1}{n^j}\right) = \sum_{n=1}^{\infty}\left(\sum_{j=2}^{\infty}\frac{(-x)^j}{j!}\frac{1}{n^j}\right) = \sum_{n=1}^{\infty}\left(e^{-x/n}-1+\frac{x}{n}\right) $$

f[x_?NumericQ] := NSum[E^(-x/n) - 1 + x/n, {n, ∞}, Method -> "EulerMaclaurin"]
g[x_] = x (Log[x] + 2 EulerGamma - 1) + 1/2;

Plot[f[x] - g[x], {x, 0, 1000}]

enter image description here

The remaining difference between the expressions is due to numerical inaccuracies and can be fixed by increasing the WorkingPrecision.

$\endgroup$
1
  • 2
    $\begingroup$ It might be a good idea to use Internal`Expm1[] here to guard against cancellation: NSum[Internal`Expm1[-x/n] + x/n, {n, 1, ∞}, Method -> "EulerMaclaurin"] $\endgroup$ Dec 13, 2021 at 9:47
0
$\begingroup$

Another way using:

$$\sum _{n=1}^{\infty } \left(\exp \left(-\frac{x}{n}\right)-1+\frac{x}{n}\right)=\int_0^{\infty } \left(\frac{x}{-1+e^z}-\frac{\sqrt{x} J_1\left(2 \sqrt{x} \sqrt{z}\right)}{\left(-1+e^z\right) \sqrt{z}}\right) \, dz$$

 f[x_] := NIntegrate[InverseLaplaceTransform[Exp[-x/n] - 1 + x/n, n, z]/(Exp[z] - 1), {z, 0, Infinity}, Method -> "LocalAdaptive"];

 fApprox[x_] := x (Log[x] + 2 EulerGamma - 1) + 1/2;
 plot1 = DiscretePlot[f[x], {x, Round[(E - 1)^Range[10, 17]]}];
 plot2 = Plot[fApprox[x], {x, 50, 10000}];
 Show[plot1, plot2, PlotRange -> All]
 DiscretePlot[f[x] - fApprox[x], {x, Round[(E - 1)^Range[10, 17]]}]
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.