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)?

  • $\begingroup$ 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. $\endgroup$
    – horchler
    Mar 11 '16 at 15:34

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

  • $\begingroup$ 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 $\endgroup$ Mar 11 '16 at 10:08
  • 2
    $\begingroup$ NIntegrate[] can handle it directly: NIntegrate[1/((z + 1)^2 + 4) 1/(Exp[-z] - 1), {z, -∞, 0, ∞}, Method -> PrincipalValue] $\endgroup$ Mar 11 '16 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 ?

  • $\begingroup$ 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. $\endgroup$ Mar 11 '16 at 12:25
  • $\begingroup$ @J.M. I totally agree with you. This was done on purpose. Cheers. $\endgroup$ Mar 11 '16 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$.

  • $\begingroup$ 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$. $\endgroup$ Mar 11 '16 at 12:32
  • 1
    $\begingroup$ 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$. $\endgroup$ Mar 11 '16 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.

  • $\begingroup$ Looks nice. Could you show some more working so we can see how you got $I$? $\endgroup$ Mar 11 '16 at 13:47
  • $\begingroup$ @MedullaOblongata: I have added some more info. $\endgroup$
    – user193810
    Mar 11 '16 at 14:20
  • $\begingroup$ 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. $\endgroup$ Mar 13 '16 at 13:18
  • $\begingroup$ 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. $\endgroup$
    – user193810
    Mar 14 '16 at 8:01
  • $\begingroup$ 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. $\endgroup$ Mar 14 '16 at 8:26

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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