I was wondering if there were additional packages to enhance Mathematica's integration capabilities - for example by including more of Gradshteyn/Ryzhik/Zwillinger tables?

Here is an example that Mathematica 10 on a desktop does not seem to be able to handle, but which can be readily found in Gradhsteyn and Ryzhik (4.267 9 or 4.267 16 in the 7th edition of Tables of Integrals, Series and Products, 2007):

Assuming[{a>0,b>0},Integrate[(t^(a - 1) - t^(b - 1))/(1 + t) 1/Log[t], {t, 0, 1}]]

Here is another example (thanks Emilio Pisanty):

Integrate[Erf[a + b*x] E^(-p*x^2), {x, -∞, ∞}]

which integrates via 8.259 1, but leaves Mathematica 10.0.2 hanging.

You should include the assumptions given in Gradhsteyn and Ryzhik (4.267 16)

\$Version

"11.1.1 for Mac OS X x86 (64-bit) (April 18, 2017)"

expr1 = Assuming[{a > 0, b > 0, r > 0},
Integrate[(t^(a - 1) - t^(b - 1))/(1 + t^r) 1/Log[t], {t, 0, 1}] //
Simplify]

(*  (1/(2*r))*
(2*r*Log[Gamma[1 + b/(2*r)]*
Gamma[(a + r)/(2*r)]] -
2*r*Log[Gamma[1 + a/(2*r)]*
Gamma[(b + r)/(2*r)]] -
a*StieltjesGamma[1,
1 + a/(2*r)] +
b*StieltjesGamma[1,
1 + b/(2*r)] +
a*StieltjesGamma[1, a/(2*r)] -
b*StieltjesGamma[1, b/(2*r)])  *)


Unfortunately, Mathematica doesn't provide the same form for an answer nor readily show that they are equivalent.

expr2 = Log[(Gamma[(a + r)/(2 r)] Gamma[b/(2 r)])/
(Gamma[(b + r)/(2 r)] Gamma[a/(2 r)])];

Assuming[{a > 0, b > 0, r > 0}, expr1 == expr2 // FullSimplify]

(*  2*r*Log[b/a] + b*StieltjesGamma[1,
1 + b/(2*r)] +
a*StieltjesGamma[1, a/(2*r)] ==
a*StieltjesGamma[1,
1 + a/(2*r)] +
b*StieltjesGamma[1, b/(2*r)]  *)


Even setting the parameters to the conditions for (4.267 9), i.e., r==1

Assuming[{a > 0, b > 0}, expr1 == expr2 /. r -> 1 // FullSimplify]

(*  2*Log[a] + a*StieltjesGamma[1,
1 + a/2] + b*StieltjesGamma[1,
b/2] == 2*Log[b] +
a*StieltjesGamma[1, a/2] +
b*StieltjesGamma[1, 1 + b/2]  *)


However, they can be shown to be numerically equivalent.

SeedRandom[0]

And @@ Table[
expr1 == expr2 /.
Thread[{a, b, r} -> RandomReal[1000, 3, WorkingPrecision -> 25]] //
Chop[#, 10^-20] &, 100]

(*  True  *)


EDIT: For r==1 (G&R 4.267 9)

expr3 = Assuming[{a > -1, b > -1},
Integrate[(t^a - t^b)/(1 + t) 1/Log[t], {t, 0, 1}] // Simplify] /. {a ->
a - 1, b -> b - 1}

(*  Log[(Gamma[(1 + a)/2] Gamma[b/2])/(Gamma[a/2] Gamma[(1 + b)/2])]  *)

expr3 == expr2 /. r -> 1

(*  True  *)

• Thanks for this great answer. It does look from other things I'm reading on the internet that Mathematica's integrator has some serious gaps. I would have expected them to include Gradhsteyn and Ryzhik as a lookup table in complement to whatever algorithm they use. – Frank Jun 6 '17 at 16:18
• Also, mine (without r) doesn't work in 10, even with the assumptions a > 0 and b > 0. Can you work verify if it works in 11.1.1? – Frank Jun 7 '17 at 0:43
• @Frank - evaluation of expr3 (edit above) works in either v11.1.1 or v10.4.1 – Bob Hanlon Jun 7 '17 at 4:16
• Interestingly, expr3 succeeds in 10.0 too, but fails if the exponents are stated as a-1, b-1. Very annoying :-) and thanks for checking! – Frank Jun 7 '17 at 4:35