# Why doesn't Mathematica solve the integral?

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


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

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.

• In my opinion, this type of question is not appropriate for this QA site. Mathematica can't solve it (not even if you put concrete values for a and b and make the integral indefinite or substitute m = Infinity). You can't solve it. Most integrals simply can't be expressed in terms of familiar functions, so this is not at all surprising. IMO any such question without an argument for why one should expect Mathematica to be able to solve a known-to-be-hard problem should be off-topic here. Jul 15 at 9:20
• @Szabolcs Thank you so much for the remark. While searching, I actually found related questions. For example, mathematica.stackexchange.com/questions/225636/… Solving this equation is very important for my research, so I am just seeking if anyone could help me solving the question.
– MK71
Jul 15 at 9:42
• Since your goal is to get a closed-form solution, not to get a closed-form solution specifically with Mathematica, math.stackexchange.com may be a better place to ask. I do not know if the question will be well-received there for the same reasons I mentioned in the comment above. When dealing with such problem, I find it useful to think about whether I truly need a closed-forum solution. Perhaps you can achieve whatever your final goal is without a closed-forum solution. Perhaps a even a numerical solution is good. Jul 15 at 9:42
• I know that there are other similar questions, but I still maintain that these are not a good fit here. I did not vote to close, as this is at the moment merely my personal opinion. As far as I am aware, there is no policy that they are off-topic (though I think there should be one). Jul 15 at 9:44
• I understand. Thank you so much for your time and suggestion.
– MK71
Jul 15 at 9:47

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.

• Thank you so much!
– MK71
Jul 15 at 19:44
• You might want to check it with some values and compare to NIntegrate just to verify it works. It might have worked in 8.0 because it was a bug and the result is incorrect - not saying this is wrong, but it seems odd that more recent versions cannot solve it. --- edit: just checked on my end and it does seem to be correct for a range of different $a$ values and matches up with NIntegrate. Jul 15 at 20:18
• It may be true that Integrate in version 8.0 has several bugs, but on the other side i remember a few cases where v.8.0 found the right solution and later versions didn't, means they had lost some knowledge. ------Since numerical integration here for several a yields the same, this analytical solution may be regarded as true. Jul 16 at 10:02
• @Akku14 just checked 12.3.1 and it's working so guess there's no regression after all. Jul 16 at 12:54

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]

• Thanks you very much for the answer!
– MK71
Jul 15 at 19:45

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.

a = 7/10; Plot3D[(1 - y*E^(b*y)*Gamma[0, y/b])/(1 + a*y), {y, 0,