Accurately computing $\sum_{j=2}^\infty \frac{(-x)^j}{j!} \zeta(j)$

Question

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]

3. 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)

Answer 1

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]

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]

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]

Answer 2

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]

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]

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.

Answer 3

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}]

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

Answer 4

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]]}]