2

我正在尝试使用 Sympy 的符号积分来找到定积分的封闭形式。特别是,我跑

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

不幸的是,在 Python 2.7.11 中,我的 Jupyter 运行、运行、运行。也许它有助于加强假设

0<H<1/2 and x>1

但我不知道该怎么做。

备注 我还使用了 Mathematica 的符号积分功能来做到这一点,它提出了一个高斯超几何函数。不幸的是,评估该函数会返回一个复数,这在评估真正的积分时并没有真正意义。因此,我希望 SymPy 可能会有所帮助。

4

1 回答 1

0

首先,Python 2.x 中的 1/2 = 0;您应该放在from __future__ import division开头或改用 Python 3。

未使用的关键假设是 x>1。实现它的一种快速方法是引入一个临时变量 y = x-1,假设它是正数,评估并替换回来。

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

如果我不导入除法,这将在几秒钟内返回一个看起来合理的结果:

(-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)

这是对你所拥有的东西的一种改进——尽管由于 1/2=0 而当然是不正确的。不幸的是,使用正确的分数,集成无法在合理的时间内完成。

我怀疑你会从 sympy 得到比从 Mathematica 更好的结果。Mathematica 结果具有复数这一事实对于困难积分来说并不罕见。这意味着必须仔细简化输出,使用复杂函数的正确分支。Mathematica 有可能自己做一些这样的简化:我建议朝着这个方向前进,也许在Mathematica.SE网站的帮助下。

于 2016-01-11T07:53:50.900 回答