Problem summing an infinite series

Question

Bug introduced in 9.0 or earlier and fixed in 11.0.0

---

Calculating this sum on Mathematica 10.3

Sum[(-1)^(r - 1)/((a^2 + r^2)r), {r, 1, Infinity}]

gives the answer

$$-\frac{1}{2a^4}+\frac{\pi^2}{12a^2}+\frac{\pi\;\text{Csch}(a\pi)}{2a^3}$$

but this is not right, and there shouldn't be any neat answer like this in terms of elementary functions. As $a\to0$, the answer should tend to $(3/4)\zeta(3)$, which doesn't have a simple form in terms of a rational multiple of $\pi^4$, as suggested by the above formula. (Typically the answer would be expressed in terms of the digamma function.)

My question is, am I doing something wrong? How should one avoid this kind of problem and be confident that the sums returned are correct?

Answer 1

Seems like a bug to me. The following is a workaround, at least on this example.

sum0 = Sum[(-1)^(r - 1)/((1. a^2 + r^2) r), {r, 1, Infinity}]

Solve::ratnz: Solve was unable to solve the system with inexact coefficients. The answer was obtained by solving a corresponding exact system and numericizing the result. >>

(*
  1/(1. + 1. a^2) 1. *
    HypergeometricPFQ[
     {1., 1.,
      (0. + 1. I) ((0. - 1. I) + (1. + 0. I) a),
      (0. - 1. I) ((0. + 1. I) + (1. + 0. I) a)},
     {2., (0. + 1. I) ((0. - 2. I) + (1. + 0. I) a),
      (0. - 1. I) ((0. + 2. I) + (1. + 0. I) a)},
     -1.]
*)

sum = sum0 // Rationalize // Simplify
(*
  (1/(4 a^2))(Log[16] + PolyGamma[0, 1/2 - (I a)/2] - 
    PolyGamma[0, 1 - (I a)/2] + PolyGamma[0, 1/2 + (I a)/2] - 
    PolyGamma[0, 1 + (I a)/2])
*)

Limit[sum, a -> 0]
% - 3/4 Zeta[3] // FullSimplify
(*
  1/16 (-PolyGamma[2, 1/2] + PolyGamma[2, 1])
  0
*)

Answer 2

Any of the methods for the Regularization option to Sum provide the same result with the expected limit.

Union[FullSimplify[
    Sum[(-1)^(r - 1)/((a^2 + r^2) r), {r, 1, Infinity}, 
     Regularization -> #]] & /@
  {"Abel", "Borel", "Cesaro", "Dirichlet", 
   "Euler"}]

(*  {(1/(4*a^2))*(HarmonicNumber[
          -(1/2) - (I*a)/2] - 
        HarmonicNumber[-((I*a)/2)] - 
        HarmonicNumber[(I*a)/2] + 
        HarmonicNumber[(1/2)*I*
            (I + a)] + Log[16])}  *)

Limit[%, a -> 0] // FunctionExpand // Simplify

(*  {(3*Zeta[3])/4}  *)

sum[a_] = (1/(4*a^2))*(
    HarmonicNumber[-(1/2) - (I*a)/2] -
     HarmonicNumber[-((I*a)/2)] -
     HarmonicNumber[(I*a)/2] +
     HarmonicNumber[(1/2)*I*(I + a)] + Log[16]);

Plot[sum[a], {a, -5, 5},
 Epilog -> {Red, AbsolutePointSize[5],
   Point[{0, 3 Zeta[3]/4}]}]

sum[a] == sum[-a] // Simplify

(*  True  *)