Mathematics Stack Exchange is a question and answer site for people studying math at any level and professionals in related fields. It's 100% free, no registration required.

Sign up
Here's how it works:
  1. Anybody can ask a question
  2. Anybody can answer
  3. The best answers are voted up and rise to the top

This is a follow-up question to my post on Stack Overflow. I want to (either analytically or numerically) integrate:

$I=\displaystyle\int_{-\infty}^{\infty}\dfrac{1}{(z+1)^2+4} \dfrac{1}{\exp(-z)-1} dz$

using MATLAB, but it tells me that the integral may not exist – the integral is undefined at $z=0$. The Cauchy principal value doesn't seem to exist either, so what does this tell us about the integral? Does it mean we can't evaluate $I$ (numerically or otherwise)?

share|cite|improve this question
Since you're using Matlab, you can compute this directly using Symbolic math and arbitrary precision: syms z; f=1/(((z+1)^2+4)*(exp(-z)-1)); I=vpa(int(f,z,-Inf,Inf)). Use double(I) to convert the result back to floating point if you like. – horchler Mar 11 at 15:34
up vote 5 down vote accepted

If you expand $\frac{1}{\exp(-z)-1}$ in a Taylor series around $z=0$, you notice that $$ \frac{1}{\exp(-z)-1}\approx -\frac{1}{z}-\frac{1}{2}-\frac{z}{12}+\ldots $$ so the function has indeed a non-integrable singularity at $z=0$. One way to assign a finite result to your integral is to subtract and re-add the singular part to your integral as $$ \int_{-\infty}^{\infty}\dfrac{1}{(z+1)^2+4} \left[\dfrac{1}{\exp(-z)-1}+\frac{1}{z}-\frac{1}{z}\right]= $$ $$ \int_{-\infty}^{\infty}\dfrac{1}{(z+1)^2+4} \left[\dfrac{1}{\exp(-z)-1}+\frac{1}{z}\right]-\mathcal{P}\int_{-\infty}^{\infty}\dfrac{1}{(z+1)^2+4} \frac{1}{z}\ , $$ where $\mathcal{P}$ is Cauchy's principal part. The final result is $-0.383448...$ (checked and agreed with Claude's answer below.)

share|cite|improve this answer
Got it in Mathematica. NIntegrate[ 1/((z + 1)^2 + 4)*(1/(Exp[-z] - 1) + 1/z), {z, -Infinity, Infinity}] - Integrate[1/((z + 1)^2 + 4)*(1/z), {z, -Infinity, Infinity}, PrincipalValue -> True] // N – Pierpaolo Vivo Mar 11 at 10:08
NIntegrate[] can handle it directly: NIntegrate[1/((z + 1)^2 + 4) 1/(Exp[-z] - 1), {z, -∞, 0, ∞}, Method -> PrincipalValue] – J. M. Mar 11 at 12:03

Noticing that the antiderivative does not exist, only numerical methods could be used.

As Ander Biguri answered your previous question, I tried to compute, as accurately as I could $$I_k=\int_{-\infty }^{-10^{-k} } \frac{dz}{\left(e^{-z}-1\right) \left((z+1)^2+4\right)}+\int_{10^{-k} }^{\infty } \frac{dz}{\left(e^{-z}-1\right) \left((z+1)^2+4\right)}$$ The results are given below $$I_1=-0.3794424331$$ $$I_2=-0.3830481054$$ $$I_3=-0.3834081111$$ $$I_4=-0.3834441111$$ $$I_5=-0.3834477111$$ $$I_6=-0.3834480711$$ $$I_7=-0.3834481071$$ $$I_8=-0.3834481107$$ $$I_9=-0.3834481110$$


After Pierpaolo Vivo's answer, it is funny to notice that what he wrote gives $$I=-0.6976073764+\frac \pi {10}$$ Who did expect $\pi$ to be here ?

share|cite|improve this answer
In principle, one could use a convergence acceleration method like Richardson extrapolation on the sequence generated in this manner, long as you pick an appropriate auxiliary sequence. – J. M. Mar 11 at 12:25
@J.M. I totally agree with you. This was done on purpose. Cheers. – Claude Leibovici Mar 11 at 19:20

Another slick way to deal with this that also exploits Pierpaolo's observation is to cancel out the singular part through an even transformation:

$$\int_{-\infty}^\infty f(x)\mathrm dx=\int_{-\infty}^\infty \frac{f(x)+f(-x)}{2}\mathrm dx$$

In this case, you should consider the integral

$$\int_{-\infty }^\infty \frac{(z+2-(z-2)\exp z)\frac{z}{\exp z-1}-5}{2\left(z^4+6z^2+25\right)} \mathrm dz$$

where the integrand has a removable singularity at $z=0$. If your quadrature routine does not sample exactly at that point, you should be good to go. (Alternatively, there are ways to handle the factor $\frac{z}{\exp z-1}$.) Mathematica in particular readily yields the result $-0.38344811106405680779$.

share|cite|improve this answer
That's a great method of solving the problem analytically. Since $z=0$ is a removable singularity, does it mean it corresponds to zero residue? i.e. the only poles we need to consider are at $z^4+6z^2+25$. – Medulla Oblongata Mar 11 at 12:32
Yes, if you perform the series expansion of that integrand at $z=0$, you should see something like $-\frac1{50}-\frac{16z^2}{1875}+\cdots$. – J. M. Mar 11 at 12:37

I did a (non-rigorous) analysis based on residuals on the plane Im[z]>0.

Let $C$ be the contour that runs from $-a$ to $a$ (in a straight line over the real axis), and back over the semi-circle $|z|=a$ with $Im[z]\ge0$, see the image (taken from Wikipedia) below.

Inside this contour, the integrand is analytical, except for a finite number of poles: the pole on $z=-1+2i$ due to the first factor, and the poles on $z=2\pi n$ for $n=1,2,...,N$ due to the second factor. Here, $N$ is the biggest integer such that $2\pi N<a$. There is also one pole exactly on the contour, at $z=0$.

According to the residue theorem, the contour-integral over C can be calculated by summing the residuals of these poles, and multiplying this by $2\pi i$. (Don't forget, here and later, that the pole at $z=0$ only counts half, because the contour goes right through it.)

Now, let $a$ go to infinity. This also means that $N$ will go to infinity. If the contribution of the semicircle part of the contour goes to zero, all that is left is $I$. (I have no proof that the contribution goes to zero. In fact, I have some doubts for $z\rightarrow i\infty$... But since the result later on seems to agree to numerical values, I ignore this, and I guess it could be made rigorously by restricting $R$ to values such that the contour falls in between the poles and not on the poles.)

The poles at $z=2n\pi i$ generate an infinite sum, that converges. I put this sum in Wolfram Alpha, and this showed that the sum can be written with digamma functions. The pole at $z=-1+2i$ gives one extra term. Altogether, the result is:

$$I=-\frac{\pi}{5} i -\frac{i}{4}\left(\psi\left(1-\frac{(1+\frac{i}{2})}{\pi}\right)-\psi\left(1+\frac{(1-\frac{i}{2})}{\pi}\right)\right)-\frac{\pi}{2-2e^{1-2i}}.$$

where $\psi$ is the digamma function. It can be simplified a little bit further by using the recurrence relation $\psi(1+x)=\psi(x)+1/x$ and the reflection formula $\psi(1-x)=\pi\cot(\pi x)+\psi(x)$.

Although the expression for $I$ has the imaginary unit $i$ in it, it is a real number.

According to wolfram alpha, it evaluates to -0.383448111064056808053521526581710702474772465557555272482. This corresponds to the other answers found numerically here, so I trust that this analysis can be made rigorously.

I tried if the digamma terms could be simplified, as there seems to be a nice symmetry in it, but I could not see how to really make it more simple.

share|cite|improve this answer
Looks nice. Could you show some more working so we can see how you got $I$? – Medulla Oblongata Mar 11 at 13:47
@MedullaOblongata: I have added some more info. – Pakk Mar 11 at 14:20
According to Mathematica, your $I=0.383448 + 1.66533*10^{-16} i$ (so it doesn't agree with the other results), and the imaginary part may not be trivial if we used different constants in my original expression. – Medulla Oblongata Mar 13 at 13:18
You are right, I missed a minus sign. But the imaginary part will always be zero if you keep on using real constants in your original expression... Because your integral will be real. – Pakk Mar 14 at 8:01
According to Mathematica, the first line of your (corrected) $I$ is $-0.383448 - 1.66533*10^{-16} i$, and the second line is $0.558761 - 5.55112*10^{-17} i$. I understand your working, but something's not quite right. – Medulla Oblongata Mar 14 at 8:26

Your Answer


By posting your answer, you agree to the privacy policy and terms of service.

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