Why doesn't Mathematica solve the integral?

Question

I am trying to solve the following integral using Mathematica,

Assuming[a > 0, Limit[Integrate[(1 - y*E^(b*y)*Gamma[0, y/b])/(1 + a*y), {y, 0, m}], m ->Infinity]]

\begin{equation}\mathop {\lim }\limits_{m \to \infty } \int_0^m {\frac{{1 - y{e^{by}}{\rm{Gamma}}\left[ {0,\frac{y}{b}} \right]}}{{1 + ay}}} dy\end{equation}

Mathematica doesn't provide any result to the integral, just shows the original equation as the output. I was trying to solve it myself, however found it too difficult to proceed further. I would be grateful if anyone could help me solving the integral.

Answer 1

You can get an analytical solution for b == 1 and a >0 with MMA version 8.0.

ii[y_, a_, b_] = (1 - y*E^(b*y)*Gamma[0, y/b])/(1 + a*y)

Integrate[ii[y, a, 1] // FunctionExpand, {y, 0, Infinity}, 
    Assumptions -> a > 0]

(*   (1/(12 a^3))E^(-1/
      a) (-12 HypergeometricPFQ[{1, 1, 1}, {2, 2, 2}, 1/a] + 
      a (-6 EulerGamma^2 + 5 \[Pi]^2 - 
  12 a E^(1/a) (EulerGamma - Log[a]) + 
  12 ExpIntegralEi[1/a] (EulerGamma - Log[a]) + 
  12 EulerGamma Log[a] - 6 Log[a]^2))   *)

For all other b you get "Integral does not converge ..."

Have no time to show it exactly, but for numbers 0<b<1 the function goes like 1/y to infinity, therefore integral is unlimited. And for b>1 it explodes to the negative.

Answer 2

Too long for a comment. Something can be done. Though the integral under consideration likely cannot be expressed in a closed form, we can consider the asymptotic of the integrand at infinity by

Normal[Series[(1 - y*E^(b*y)*Gamma[0, y/b])/(1 + a*y), {y, Infinity, 
4}, Assumptions -> a > 0 && b > 0]] // Simplify

(1/(a^4 y^4))(-1 + a y - a^2 y^2 + a^3 y^3 + b E^(-(y/b) + b y) (1 + a (b - y) + a^2 (2 b^2 - b y + y^2) + a^3 (6 b^3 - 2 b^2 y + b y^2 - y^3)))

and integrate it

Integrate[%, {y, 10, m}]

ConditionalExpression[ 1/a^4 (1/3000 (-1 - 300 a^2 - 5 E^(-(10/b) + 10 b) + ( 50 E^(-(10/b) + 10 b))/ b - (99 + 20 a - 550 a^2) b E^(-(10/b) + 10 b) + (5 - 249 a - 25 a^2 + 900 a^3) b^2 E^(-(10/b) + 10 b) - (-50 - 5 a + 348 a^2 + 60 a^3) b^3 E^(-(10/b) + 10 b) - 2 a (-25 - 5 a + 447 a^2) b^4 E^(-(10/b) + 10 b) + 10 a^2 (10 + 3 a) b^5 E^(-(10/b) + 10 b) + 300 a^3 b^6 E^(-(10/b) + 10 b) + 5 a (3 + 40 E^(-(10/b) + 10 b)) - 1/b^2 500 (a b (-4 + b^2) (-1 + b^2)^2 + (-1 + b^2)^3 + 6 a^3 b^3 (-4 + 6 b^2 - 4 b^4 + b^6) + a^2 b^2 (-11 + 18 b^2 - 9 b^4 + 2 b^6)) ExpIntegralEi[( 10 (-1 + b^2))/b] - 3000 a^3 Log[10]) + 1/(6 b^2 m^3) ((a b (-4 + b^2) (-1 + b^2)^2 + (-1 + b^2)^3 + 6 a^3 b^3 (-4 + 6 b^2 - 4 b^4 + b^6) + a^2 b^2 (-11 + 18 b^2 - 9 b^4 + 2 b^6)) m^3 ExpIntegralEi[((-1 + b^2) m)/b] - b (E^(((-1 + b^2) m)/b) m^2 + 6 a^3 b^7 E^(((-1 + b^2) m)/b) m^2 + 2 a^2 b^6 E^(((-1 + b^2) m)/b) m (3 a + m) + b^2 E^(((-1 + b^2) m)/b) (2 - 4 a m - 2 m^2 + 11 a^2 m^2) + b (-2 - E^(((-1 + b^2) m)/b) m - 6 a^2 m^2 + a m (3 + 4 E^(((-1 + b^2) m)/b) m)) + b^3 E^(((-1 + b^2) m)/ b) (m - 5 a^2 m + 18 a^3 m^2 + a (2 - 5 m^2)) - a b^5 E^(((-1 + b^2) m)/ b) (-2 a m - m^2 + 6 a^2 (-2 + 3 m^2)) - b^4 E^(((-1 + b^2) m)/ b) (-a m + 12 a^3 m - m^2 + a^2 (-4 + 7 m^2)) - 6 a^3 b m^3 Log[m]))), Re[m] > 0 || m \[NotElement] Reals]

Now

Limit[%, m -> Infinity]

ConditionalExpression[-\[Infinity], a b > 0 && b^3 > b && E^(1/b) != 0]

produces a certain condition on the parameters when the limit is infinite. Taking {y, Infinity, 5} in the above, a very close result is obtained

ConditionalExpression[-\[Infinity], b > 1 && E^(1/b + 2 b^3) != 0 && a b^2 E^(10 b^3) > 0]

Answer 3

The integral appears to be solvable in Mathematica 8.0.1 when substituting a and b with b=Prime[3] and a=Prime[5] and applying integration by parts 3 times:

Part1 = (1 - y*E^(Prime[3]*y)*Gamma[0, y/Prime[3]]);
Part2 = 1/(1 + Prime[5]*y);
u = Part1;
du = D[Part1, y];
v = Integrate[Part2, y];
dv = Part2;
u*v - Integrate[v*du, y]
Integrate[Part1*Part2, y];

which gives:

(1/11)*Log[11*y + 1]*(1 - E^(5*y)*y*Gamma[0, y/5]) - 
   (1/11)*Integrate[Log[11*y + 1]*(E^((24*y)/5) - 
            E^(5*y)*Gamma[0, y/5] - 5*E^(5*y)*y*Gamma[0, y/5]), y]

Integration by parts of the unsolved expression within the integral:

Part2 = (E^(24 y/5) - E^(5 y) Gamma[0, y/5] - 
    5 E^(5 y) y Gamma[0, y/5]) ;
Part1 = Log[1 + 11 y];
u = Part1;
du = D[Part1, y];
v = Integrate[Part2, y];
dv = Part2;
u*v - Integrate[v*du, y]

gives:

11*Integrate[(E^(5*y)*y*Gamma[0, y/5])/(1 + 11*y), y] - 
   E^(5*y)*y*Gamma[0, y/5]*Log[1 + 11*y]

Integrating by parts of the unsolved expression a last time:

Part2 = E^(5 y) y Gamma[0, y/5];
Part1 = 1 + 11 y;
u = Part1;
du = D[Part1, y];
v = Integrate[Part2, y];
dv = Part2;
u*v - Integrate[v*du, y]

gives:

(-(11/600))*((365/24)*E^((24*y)/5) - (48/5)*ExpIntegralEi[(24*y)/5] - 
        24*y*ExpIntegralEi[(24*y)/5] + (24/5)*E^(5*y)*(-2 + 5*y)*
     Gamma[0, y/5]) + 
   (1 + 11*y)*((1/24)*E^((24*y)/5) - (1/25)*ExpIntegralEi[(24*y)/5] + 
        (1/25)*E^(5*y)*(-1 + 5*y)*Gamma[0, y/5])

which has no unsolved integral expressions in it, so therefore it is possible that your integral is solvable, but you need to refine this answer by including not only the unsolved expressions and also try different values of a and b to make it exact. Also since Mathematica 8 is an old version, newer versions might give different answers. That because of previous programming bugs in integration in Mathematica 8.

Answer 4

A plot of the integrand shows the region where it gets very (infinitely) large so the general integral will have a hard time...:

a = 7/10; Plot3D[(1 - y*E^(b*y)*Gamma[0, y/b])/(1 + a*y), {y, 0, 
12000}, {b, 0.1, 1.5}, AxesLabel -> Automatic, 
WorkingPrecision -> 30, PlotPoints -> 500, PlotRange -> {-10^5, 10}]