I have to finish following integral $$\int_0^{+\infty} \log_2(1+a x) (1-e^{-\frac{x}{2}})^{b-1} \frac{1}{2}e^{-\frac{x}{2}} dx$$ After a week of work on this, I found it may be impossible to have a integral analytically. So I think I need to have an approximation of this.

The shape of the integrated function is as following. With fixed b and varying a: With fixed b and varying a With fixed a and varying b:With fixed a and varying b

Could you give me a hint on the approximation of the integral?


Solution 1:

You can approximate the integral using Gauss-Laguerre quadrature (http://en.wikipedia.org/wiki/Gauss%E2%80%93Laguerre_quadrature) for integrals of the form $$ \int_0^{\infty} f(x)e^{-x}dx \approx \sum_{i=1}^n w_i f(x_i), $$ where $x_i,i=1,...,n$ is the root of the $n$-th degree Laguerre polinomial $L_n(x)$ and $ w_i,i=1,...,n$ is $$ w_i=\frac{x_i}{(n+1)^2(L_{n+1}(x_i))^2}. $$ Picking $n=1$ results in an approximation which, albeit crude, has the advantage of being a simple (and interpretable) function of $a$ and $b$, as $L_1(x)=-x+1,L_2(x)=1/2(x^2-4x+2),x_1=1$ and $w_1=1$. Defining $f(u)=\log_2(1+2au)(1-e^{-u})^{b-1}$, the integral can be approximated as follows: $$ \int_0^{\infty} f(u)e^{-u}du\approx 1\cdot f(1)=\log_2(1+2a)(1-e^{-1})^{(b-1)}\quad (*) $$ This table shows a few instances of the correct integral $I(a,b)$ and approximation $\hat I(a,b)$.

a b I I^
0 -1 0 0 
0 0 0 0 
0 1 0 0 
1 0 2.650694 2.507374 
1 1 1.331479 1.584963 
1 2 0.9013049 1.001887 
2 1 1.934489 2.321928 
2 2 1.268749 1.467738 
2 3 0.9494982 0.9277877 
3 2 1.508186 1.774587 
3 3 1.117553 1.121753 
3 4 0.8908656 0.709083 

More accurate numerical results can be obtained with $n=2$ or larger but the simplicity of $(*)$ will be lost. Hope this helps!

Solution 2:

I know this is a sort of cold-case answer but I have just found a different method and I think the idea may be of some interest.

Intro: it is "easy" to numerically compute the integral for given values of $a$ and $b$, just type

f <- function(x,a=1,b=1) log(1+a*x,2)*(1-exp(-x/2))**(b-1)*exp(-x/2)/2
integrate(f,0,Inf,a=2,b=3)

to get "0.9494982 with absolute error < 9.6e-07", say when $a=2,b=3$, using R (many other programs would, of course, do the job). However, the question is asking a simple approximation. Assuming this means a "simple formula" you can try a different approach.

Eureqa: generate a grid of $\left(a_i,b_i,\int_0^{\infty} f(x,a_i,b_i)dx\right)$ and let this symbolic regression software find the approximate formula for you. The idea is taken from this paper by Stoutemyer: I do not want to spoil the fun but, briefly, the paper explains how this software can approximate relationships between variables, guessing the unknown functional form and can even find the exact formula when there is one, sometimes improving the results of symbolic Computer Algebra Systems. The software can be downloaded here

Hence, I computed the value of the integral on an evenly spaced grid of $11\times11$ couples $(1\leq a_i\leq 10,1\leq b_i\leq 10)$. Then I imported these 441 values in Eureqa and let it work for about 5 minutes to obtain the approximation $$ \frac{11.21 + 6.537ab}{9.084b + b^2 + ab^2}. $$ This is simple and accurate (in the previously described domain of $a$ and $b$, the absolute error is often smaller than 0.01, with larger errors clustered maybe where $a,b\approx 1$, see Figure). Perhaps this is not the place to discuss other more complex approximations involving logs provided by the software, different domain of interest or optimal Padè (rational) approximations: all of this can be done in Eureqa but ultimately depends on what the results are needed for.

But I think, as pointed out by Stoutemyer, that the idea to use a numerical integrator to generate data to be used to find approximate or exact relationships by symbolic regression is intriguing as, among other things, the software is picking the functional form (largely) by itself.

Absolute error for $0\leq a,b\leq 10$