5

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

  • 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 '16 at 15:34

4 Answers4

5

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

  • 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 '16 at 10:08
  • 2
    NIntegrate[] can handle it directly: NIntegrate[1/((z + 1)^2 + 4) 1/(Exp[-z] - 1), {z, -∞, 0, ∞}, Method -> PrincipalValue] – J. M. ain't a mathematician Mar 11 '16 at 12:03
5

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$$

Edit

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 ?

4

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

  • 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 '16 at 12:32
  • 1
    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. ain't a mathematician Mar 11 '16 at 12:37
4

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.

  • Looks nice. Could you show some more working so we can see how you got $I$? – Medulla Oblongata Mar 11 '16 at 13:47
  • @MedullaOblongata: I have added some more info. –  Mar 11 '16 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 '16 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. –  Mar 14 '16 at 08:01
  • According to Mathematica, the first line of your (corrected) $I$ is $-0.383448 - 1.6653310^{-16} i$, and the second line is $0.558761 - 5.5511210^{-17} i$. I understand your working, but something's not quite right. – Medulla Oblongata Mar 14 '16 at 08:26
  • Something is wrong indeed. I don't have the time to find out what exactly (probably it was just one minus sign), so I took the second line out and mentioned the steps that could be taken to simplify further. By the way: try calculating at higher precision in Mathematica. The imaginary part of your approximation is in the order of the machine epsilon... –  Mar 14 '16 at 09:40