Erroneous Wolfram result for $\sum_{k=1}^\infty (k^3 + a^3)^{-1}$, looking for correct formula I was trying to get some interesting result for $\zeta(3)$, exploring the following function:
$$W(a) = \sum_{k=1}^\infty \frac{1}{k^3 + a^3}, \mbox{ with } \lim_{a\rightarrow 0} W(a) = \zeta(3).$$
Let $w_1, w_2, w_3$ be the three roots (one real, two complex) of $(w+1)^3+a^3=0$, with $w_1=-(a+1)$. Also, $a$ is a real number. Using Wolfram Alpha (see here), I get
$$W(a)=\frac{-1}{3}\cdot\sum_{j=1}^3 W_j(a), \mbox{ with } W_j(a) = \frac{\psi^{(0)}(-w_j)}{(w_j+1)^2}.$$
Here $\psi^{(0)}$ is the digamma function. The result is wrong because $W_1(a) \rightarrow \infty$ as $a\rightarrow 0^+$ while $W_2(a)$ and $W_3(a)$ remain bounded. Indeed using $a=0.0001$, Wolfram yields $W(a)\approx -2334.16$, see here. Surprising, with $a=0.01$ it yields $W(a)\approx 1.20206$ which is very close to the true result.
Surprisingly, Wolfram knows (see here) that
$$\lim_{a\rightarrow 0} W(a) = -\frac{\psi^{(2)}(1)}{2}.$$
Of course (this is a well known fact), $\zeta(3)=-\psi^{(2)}(1)/2$ and thus Wolfram is correct this time.
My question:
What is going on with this computation (or is it me?), and what is the correct formula for $W(a)$?
Update
See the two answers below proving that I was wrong, and that the Mathematica formula I though was incorrect, is indeed right. Kudos Mathematica! You were successful at solving a nice problem involving a few challenging steps, and coming with a somewhat unexpected but neat formula involving derivatives of the digamma function instead of the classic $\zeta(3)$.
Final note
It is possible to use a different, simpler approach that does not involve complex numbers. Consider
$$V(a) =\sum_{k=1}^\infty (-1)^{k+1}\frac{1}{k(k^2-a^2)}.$$
Wolfram is able to compute the limit of $V(a)$ as $a\rightarrow 0$, and returns the correct value $3\zeta(3)/4$, see here. It is easy to establish that
$$V(a)=\frac{1}{a^2} \Big[\int_0^\infty \frac{\cosh(ax)}{1+e^x} dx -\log 2\Big].$$
To compute $\lim_{a\rightarrow 0} V(a)$, we apply L'Hospital Rule twice to the above expression, the denominator in this case being $a^2$. This yields
$$\lim_{a\rightarrow 0}V(a) = \frac{1}{2}\lim_{a\rightarrow 0}\int_0^\infty \frac{x^2\cosh(ax)}{1+e^x}dx =\frac{1}{2}\int_0^\infty \frac{x^2}{1+e^x}dx=\frac{3\zeta(3)}{4}.$$
Here, we assume $a<1$.
 A: We have the partial fraction decomposition
$$\frac{ca^2}{k^3+a^3}=\frac{-\omega}{k-a/\omega }+\frac{\omega -1}{k+a}+\frac{1}{k-a \omega},$$
where $c:=3(\omega-1)$ and $\omega:=e^{i\pi/3}$.
Also,
$$\sum_{k=1}^n\frac1{k+b}=\ln n-\psi(1+b)+o(1)$$
(as $n\to\infty$), where $\psi$ is the digamma function.
Collecting the pieces, for $a\in(-1,\infty)\setminus\{0\}$ we get
$$s(a):=\sum_{k=1}^\infty\frac1{k^3+a^3}
=\frac1{ca^2}\,
\left((1-\omega) \psi(1+a)+\omega\psi\left(1-a/\omega\right)
-\psi(1-a \omega)\right).$$
For $a\to0$,
$$s(a)=-\frac{\psi ^{(2)}(1)}{2}-\frac{\pi ^6 a^3}{945}+O\left(a^4\right)
=\zeta(3)-\frac{\pi ^6 a^3}{945}+O\left(a^4\right).$$
Here is the graph $\{(a,s(a))\colon0<a\le1\}$, with $s(0)=\zeta(3)=1.2020\ldots$:

(I am not geting instability.)
A: I think the statement in the OP that $W_2(a)$ and $W_3(a)$ remain bounded when $a\rightarrow 0$ is mistaken, so that there is no inconsistency with the Mathematica result.
The three roots of $(w+1)^3+a^3=0$ are
$$w_1= -a-1,\;\; w_2= \tfrac{1}{2} \left(-i \sqrt{3} a+a-2\right),\;\;w_3= \tfrac{1}{2} \left(i \sqrt{3} a+a-2\right).$$
Then the denominator $(w+1)^2$ vanishes for all three roots when $a\rightarrow 0$, while the numerator remains finite (equal to $-\gamma_{\rm Euler}$).
And indeed, a numerical check suggets that the Mathematica output is actually correct, and the erroneous numerical result for small $a$ is a numerical instability in the computation of the digamma function. See these two plots that compare the digamma expression (blue) with a numerical evaluation of the sum (gold), as a function of $a$. For $a\gtrsim 0.01$ the two answers are nearly indistinguishable.


