5

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.

Frank
  • 387
  • 1
  • 7

1 Answers1

6

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  *)
Bob Hanlon
  • 157,611
  • 7
  • 77
  • 198
  • 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 06 '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 07 '17 at 00:43
  • @Frank - evaluation of expr3 (edit above) works in either v11.1.1 or v10.4.1 – Bob Hanlon Jun 07 '17 at 04: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 07 '17 at 04:35