2

I am trying to use Sympy's symbolic integration to find a closed form for a definite integral. In particular, I run

from sympy import *
x, s, H = symbols('x s H', real=True, nonnegative=True)
integrate(1/((1-s)**(1/2-H)*(x-s)**(1/2-H)),(s,0,1))

Unfortunately, with Python 2.7.11 my Jupyter runs and runs and runs. Maybe it helps to strengthen the assumptions adding

0<H<1/2 and x>1

but I didn't find out how to do it.

Remark I have also used Mathematica's symbolic integration capabilities to do it and it has come up with a Gauss hypergeometric function. Unfortunately, evaluating that function returns a complex number which doesn't really make sense in evaluating a real integral. Hence my hope that SymPy might help.

Ben
  • 465
  • 2
  • 6
  • 17
  • 1
    Which error do you receive? Could you add the the equation to your question, just to make sure that you did not mess up something regarding the parentheses?! How do you incorporate your additional assumptions (H between 0 and 0.5 and x >1)? – Cleb Dec 15 '15 at 12:08

1 Answers1

0

First of all, 1/2 = 0 in Python 2.x; you should put from __future__ import division at the beginning or use Python 3 instead.

The key unused assumption is that x>1. A quick way to implement it is to introduce a temporary variable y = x-1, assume it to be positive, evaluate and substitute back. 

from __future__ import division
from sympy import *
s, x = symbols('s x')
y, H = symbols('x H', positive=True)
f = 1/((1-s)**(1/2-H)*(x-s)**(1/2-H))         % function to be integrated
integrate(f.subs(x,y+1), (s,0,1)).subs(y,x-1)  % sub, integrate, sub back

If I didn't import division, this would return within several seconds with a reasonable-looking result:

(-1)**H*(x - 1)**H*exp(I*pi*H)*gamma(H + 1)*hyper((-H, H + 1), (H + 2,), -exp_polar(2*I*pi)/(x - 1))/gamma(H + 2)

which is kind of an improvement on what you had -- though of course incorrect due to 1/2=0. Unfortunately, with the correct fraction in place integration fails to finish in reasonable time.

I doubt that you will get a better result from sympy than from Mathematica. The fact that the Mathematica result has complex numbers is not unusual for difficult integrals. It means one has to carefully simplify the output, using the correct branches of the complex functions. It's possible Mathematica can do some of this simplification on its own: I suggest proceeding in this direction, perhaps with the help of Mathematica.SE site.

Community
  • 1
  • 1