0

I've been trying to decipher some old (undocumented, poorly written) code and have reached a sticking point. It boils down to evaluating an integral of the form $$I = \int_{0}^{\infty}L(x)\cdot\frac{1}{1+\exp(a+bx)}dx$$ where $L(x) = \frac{1}{1+(x-c)^2}$.

The original writers then (I think) split the integral at $x_*=-\frac{b}{a}$ and apply a binomial expansion (ensuring the exponent in the denominator is less than 1) to $(1+\exp(a+bx))^{-1}$ yielding: $$I = \int_{0}^{x_*}L(x)dx +\sum_{k=1}^{\infty}(-1)^k\Big[\int_{0}^{x_*}L(x)e^{k(a+bx)}dx + \int_{x_*}^{\infty}L(x)e^{-k(a+bx)}dx\Big]$$

The bit I am stuck at is what they do next with $L(x)$. The result has some fixed constants which I can only assume are constants of integration. Taking the first part of the above expansion where $k=0$, the code seems to be doing something like the following:

$$\int_0^{x_*}L(x)dx = \sum_j\alpha_j\int_0^{x_*}\frac{e^{-\beta_j(x-c)}}{x}dx$$ for some pre-defined/computed j-dependent $\alpha, \beta$ which are a bit of a mystery to me. It looks like a bit like a series expansion making some use of the exponential integral $Ei(x)$ ?! I think this part of the integral is only not performed analytically so it's in the same format as the $k\neq0$ parts of the expansion.

My questions are:

(1) Where might these constants have come from and/or what would the next step be to compute this integral? Any theories are welcome!

(2) Could I simply abandon all these expansions and compute the integral analytically using a Laplace transform, particularly making use of the property $\int_0^{\infty}f(x)g(x)dx = \int_0^{\infty}\mathcal{L}\{f(x)\}(s)\cdot \mathcal{L}^{-1}\{g(x)\}(s)ds$ (ref: https://en.wikipedia.org/wiki/Laplace_transform#Evaluating_integrals_over_the_positive_real_axis) ?

Ideally I'd like to extend whatever evaluation method I settle on to use some different forms of $L(x)$.

Thanks in advance! :)

ahoffus
  • 15
  • 4
  • I don't know why you would do anything special to compute $\int_0^{x_*} L(x) dx$ when it is just the reciprocal of a quadratic polynomial, so you can just integrate it exactly. – Ian Jul 26 '21 at 15:59
  • @Ian Me neither! It looks like the original authors of the code have employed the same method for the $\int_0^{x_}L(x)dx$ as they have for all the $\int_0^{x_}L(x)\exp(k(a+bx))$ integrals, just with $k=0$ instead. – ahoffus Jul 26 '21 at 16:01
  • Well, when you put it that way, perhaps the point is that you can analytically determine the partial fraction decomposition of $L(x)$ which reduces all these integrals to the form you wrote. So if you have some kind of fast Ei implementation then this might be a good way to go. It still seems completely unnecessary to numerically evaluate the part with no exponential, though. – Ian Jul 26 '21 at 16:05
  • That sounds like sensible next step - I'll give that a go, thanks! Yes there's a bit of re-writing needed to avoid the unnecessary work that can be done exactly – ahoffus Jul 26 '21 at 16:09
  • Personally if I was told "evaluate this class of integrals numerically" without any relatively sophisticated software, I would write $I=b^{-1} \int_0^\infty e^{-x} L(x/b) \frac{1}{e^{-x}+e^a} dx$ and then use Gauss-Laguerre quadrature. That assumes $b>0$ of course, if $b<0$ then you have a rather different situation. – Ian Jul 26 '21 at 16:10
  • ...Hold on, can I get a sign check on the parameters $a,b,c$? What I just said was assuming they are all positive, but you seem to have $x^*>0$ which suggests exactly one of $a$ and $b$ is negative. – Ian Jul 26 '21 at 16:21
  • $c$ is positive. $a$ and $b$ could be any combination of sign. I should have written $x_* = \max(0, -b/a)$. – ahoffus Jul 26 '21 at 16:33
  • OK yeah, what I said is nonsense when $b<0$, and won't perform super great when $a<0<b$ either. In view of that, maybe this partial fractions method is better. – Ian Jul 26 '21 at 16:36

1 Answers1

1

It seems that the approach uses the fact that

$$L(x)=\frac{1}{2i} \left ( \frac{1}{x-(c+i)}-\frac{1}{x-(c-i)} \right ).$$

Next the point is that you can do a geometric series expansion of the other term, but which one is appropriate depends on which side of $x_*$ you're on:

  • When $a+bx>0$, the exponential is bigger, and so you multiply top and bottom by $\exp(-a-bx)$ to get $\frac{\exp(-a-bx)}{1+\exp(-a-bx)}$ leaving you with $\sum_{n=1}^\infty (-1)^{n-1} \exp(-n(a+bx))$
  • When $a+bx<0$, the exponential is smaller, so you just have $\sum_{n=0}^\infty (-1)^n \exp(n(a+bx))$.

Which way is on which side of $x_*$ depends on the sign of $b$, contrary to what you wrote which requires $b>0$. And of course this only even matters if $a$ and $b$ have different signs. If they have the same sign then $x_*$ is out of the domain of integration and you just do the first thing throughout.

So overall you have some case work.

If $a$ and $b$ are both positive, then

$$f(x)=\sum_{n=1}^\infty (-1)^{n-1}\frac{1}{2i} \left ( \frac{1}{x-(c+i)}-\frac{1}{x-(c-i)} \right ) \exp(-n(a+bx)).$$

If $a$ and $b$ are both negative, then

$$f(x)=\sum_{n=0}^\infty (-1)^n \frac{1}{2i} \left ( \frac{1}{x-(c+i)}-\frac{1}{x-(c-i)} \right ) \exp(n(a+bx)).$$

Otherwise:

$$f(x)=\begin{cases} \sum_{n=1}^\infty (-1)^{n-1}\frac{1}{2i} \left ( \frac{1}{x-(c+i)}-\frac{1}{x-(c-i)} \right ) \exp(-n(a+bx)) & a+bx>0 \\ \sum_{n=0}^\infty (-1)^n \frac{1}{2i} \left ( \frac{1}{x-(c+i)}-\frac{1}{x-(c-i)} \right ) \exp(n(a+bx)) & a+bx<0 \end{cases}.$$

In any case the last step to get these into the form you had up there is to simply change variables to $x-(c\pm i)$ depending on which term we're talking about.

Ian
  • 104,572