A well know result in theoretical physics is that a sum over Matsubara fermionic frequencies, i.e.:
$$ S = \sum_{n=-\infty}^{\infty} h(\omega_n) \hspace{32pt} \omega_n=(2n+1)\frac{\pi}{\beta} $$
can be rewritten as a contour integral as
$$ S = \frac{-1}{2 \pi \mathrm{i}} \oint_C \mathrm{d} z g(z) h(-\mathrm{i}z) \hspace{32pt} g(z)=\frac{\beta}{2} \tanh(\frac{\beta z}{2})$$
where $C$ is an infinite contour enclosing counterclockwise the imaginary axis in the complex plane. This result is demonstrated, for instance, in Altand & Simons, "Condensed Matter Field Theory", Cambridge University Press (2010).
If the integrand decays fast enough the contour $C$ can be split in two contours, parallel to the imaginary axis. Enough with the physics! To make a long story short the theory implies that the following functions in Mathematica should give the same result (for any choice of dist):
h[omega_] := omega^2/(omega^4 + omega^3 + 1)
omega[n_] := (2*n + 1)*\[Pi]/b
s1[bval_] := Sum[h[omega[n]], {n, -Infinity, Infinity}] /. {b -> bval}
g[z_, b_] := b/2 Tanh[z b/2]
s2[b_, dist_] := 1/(2*\[Pi])*NIntegrate[(g[z1, b]*h[-I*z1] - g[z2, b]*h[-I*z2])
/. {z1 -> I*x + dist, z2 -> I*x - dist}, {x, -Infinity, Infinity},
Method -> "DoubleExponentialOscillatory", MaxRecursion -> 16]
i.e. s1[b] should evaluate to the same result as s2[b,d] for any choice of b and d. However, while the sum in s1 is evaluated instantly, the integral in s2 fails to converge, even if a plot of the integrand looks good enough:
Note: The peaks are in infinite number, periodically spaced, but their height is rapidly decreasing. Probably this confuses Mathematica integration strategies.
I tried to change the integration method, AccuracyGoal and Precision goal, as long as increasing MaxRecursion. Sometimes the results of the integral s2 are very near to the exact result from s1, but I always get the NIntegrate::ncvb error (and most often also NIntegrate::slwcon).
I also tried to convert the infinite integral to a finite one, using:
Leading to:
s3[b_, dist_] := 1/(2*\[Pi])*NIntegrate[(g[z1, b]*h[-I*z1] - g[z2, b]*h[-I*z2])
*(1 + t^2)/(1 - t^2)^2 /. {z1 -> I*t/(1 - t^2) + dist,z2 -> I*t/(1 - t^2) - dist},
{t, -1, 1}]
And the peaks become now an oscillating behaviour at the borders, which apparently Mathematica cannot handle:
This conversion seems to help somehow, but I never get an errorless, precise answer from Mathematica.
g
is undefined in the question. I presume it should beg[z_, bval_] := bval/2 Tanh[z bval/2]
$\endgroup$