对一组特定参数使用 SymPy 符号积分时出错

Error using SymPy symbolic integration for a specific set of parameters

我使用符号积分求出以下一系列(类似的)积分的数值,因为这些函数收敛速度较慢,而且迄今为止的数值方法得出的值与符号计算的值大不相同(无论我试过什么相对和绝对容忍度):

from sympy import *
import numpy as np

a1 = 1001
a2 = 1001
a3 = 951

Delta = lambda s: ((a1**2+s)**Rational(1, 2))*((a2**2+s)**Rational(1, 2))*((a3**2+s)**Rational(1, 2))

s = Symbol('s')

I1 = re(N(2*np.pi*a1*a2*a3*integrate(1/((a1**2 + s) * Delta(s)), (s, 0, oo))))
I2 = re(N(2*np.pi*a1*a2*a3*integrate(1/((a2**2 + s) * Delta(s)), (s, 0, oo))))
I3 = re(N(2*np.pi*a1*a2*a3*integrate(1/((a3**2 + s) * Delta(s)), (s, 0, oo))))

积分适用于 a1、a2、a3 的大多数值。但是,偶然地,我发现了一组导致错误的值:a1 = 1001,a2 = 1001,a3 = 951。出现以下错误:

python3 test.py
Traceback (most recent call last):
  File "/opt/anaconda3/lib/python3.7/site-packages/sympy/core/assumptions.py", line 454, in getit
    return self._assumptions[fact]
KeyError: 'zero'

During handling of the above exception, another exception occurred:

Traceback (most recent call last):
  File "/opt/anaconda3/lib/python3.7/site-packages/sympy/core/assumptions.py", line 454, in getit
    return self._assumptions[fact]
KeyError: 'extended_real'

During handling of the above exception, another exception occurred:

Traceback (most recent call last):
  File "/opt/anaconda3/lib/python3.7/site-packages/sympy/core/assumptions.py", line 454, in getit
    return self._assumptions[fact]
KeyError: 'finite'

During handling of the above exception, another exception occurred:

Traceback (most recent call last):
  File "test.py", line 14, in <module>
    I1 = re(N(2*np.pi*a1*a2*a3*integrate(1/((a1**2 + s) * Delta(s)), (s, 0, oo))))
  File "/opt/anaconda3/lib/python3.7/site-packages/sympy/integrals/integrals.py", line 1571, in integrate
    return integral.doit(**doit_flags)
  File "/opt/anaconda3/lib/python3.7/site-packages/sympy/integrals/integrals.py", line 581, in doit
    ret = try_meijerg(function, xab)
  File "/opt/anaconda3/lib/python3.7/site-packages/sympy/integrals/integrals.py", line 553, in try_meijerg
    res = meijerint_definite(function, x, a, b)
  File "/opt/anaconda3/lib/python3.7/site-packages/sympy/integrals/meijerint.py", line 1860, in meijerint_definite
    res = _meijerint_definite_2(f, x)
  File "/opt/anaconda3/lib/python3.7/site-packages/sympy/integrals/meijerint.py", line 1969, in _meijerint_definite_2
    res = _meijerint_definite_3(g, x)
  File "/opt/anaconda3/lib/python3.7/site-packages/sympy/integrals/meijerint.py", line 1981, in _meijerint_definite_3
    res = _meijerint_definite_4(f, x)
  File "/opt/anaconda3/lib/python3.7/site-packages/sympy/integrals/meijerint.py", line 2060, in _meijerint_definite_4
    res += C*_int0oo(f1_, f2_, x)
  File "/opt/anaconda3/lib/python3.7/site-packages/sympy/integrals/meijerint.py", line 1305, in _int0oo
    return meijerg(a1, a2, b1, b2, omega/eta)/eta
  File "/opt/anaconda3/lib/python3.7/site-packages/sympy/core/decorators.py", line 266, in _func
    return func(self, other)
  File "/opt/anaconda3/lib/python3.7/site-packages/sympy/core/decorators.py", line 136, in binary_op_wrapper
    return func(self, other)
  File "/opt/anaconda3/lib/python3.7/site-packages/sympy/core/expr.py", line 266, in __truediv__
    return Mul(self, denom)
  File "/opt/anaconda3/lib/python3.7/site-packages/sympy/core/cache.py", line 72, in wrapper
    retval = cfunc(*args, **kwargs)
  File "/opt/anaconda3/lib/python3.7/site-packages/sympy/core/operations.py", line 85, in __new__
    c_part, nc_part, order_symbols = cls.flatten(args)
  File "/opt/anaconda3/lib/python3.7/site-packages/sympy/core/mul.py", line 266, in flatten
    if not a.is_zero and a.is_Rational:
  File "/opt/anaconda3/lib/python3.7/site-packages/sympy/core/assumptions.py", line 458, in getit
    return _ask(fact, self)
  File "/opt/anaconda3/lib/python3.7/site-packages/sympy/core/assumptions.py", line 513, in _ask
    _ask(pk, obj)
  File "/opt/anaconda3/lib/python3.7/site-packages/sympy/core/assumptions.py", line 501, in _ask
    a = evaluate(obj)
  File "/opt/anaconda3/lib/python3.7/site-packages/sympy/core/expr.py", line 912, in _eval_is_extended_negative
    return self._eval_is_extended_positive_negative(positive=False)
  File "/opt/anaconda3/lib/python3.7/site-packages/sympy/core/expr.py", line 870, in _eval_is_extended_positive_negative
    if self.is_extended_real is False:
  File "/opt/anaconda3/lib/python3.7/site-packages/sympy/core/assumptions.py", line 458, in getit
    return _ask(fact, self)
  File "/opt/anaconda3/lib/python3.7/site-packages/sympy/core/assumptions.py", line 513, in _ask
    _ask(pk, obj)
  File "/opt/anaconda3/lib/python3.7/site-packages/sympy/core/assumptions.py", line 501, in _ask
    a = evaluate(obj)
  File "/opt/anaconda3/lib/python3.7/site-packages/sympy/core/expr.py", line 846, in _eval_is_positive
    finite = self.is_finite
  File "/opt/anaconda3/lib/python3.7/site-packages/sympy/core/assumptions.py", line 458, in getit
    return _ask(fact, self)
  File "/opt/anaconda3/lib/python3.7/site-packages/sympy/core/assumptions.py", line 513, in _ask
    _ask(pk, obj)
  File "/opt/anaconda3/lib/python3.7/site-packages/sympy/core/assumptions.py", line 513, in _ask
    _ask(pk, obj)
  File "/opt/anaconda3/lib/python3.7/site-packages/sympy/core/assumptions.py", line 501, in _ask
    a = evaluate(obj)
  File "/opt/anaconda3/lib/python3.7/site-packages/sympy/core/expr.py", line 909, in _eval_is_extended_positive
    return self._eval_is_extended_positive_negative(positive=True)
  File "/opt/anaconda3/lib/python3.7/site-packages/sympy/core/expr.py", line 875, in _eval_is_extended_positive_negative
    n2 = self._eval_evalf(2)
  File "/opt/anaconda3/lib/python3.7/site-packages/sympy/functions/special/hyper.py", line 690, in _eval_evalf
    v = mpmath.meijerg(ap, bq, z, r)
  File "/opt/anaconda3/lib/python3.7/site-packages/mpmath/functions/hypergeometric.py", line 1058, in meijerg
    return ctx.hypercomb(h, a+b, **kwargs)
  File "/opt/anaconda3/lib/python3.7/site-packages/mpmath/functions/hypergeometric.py", line 127, in hypercomb
    [ctx.rgamma(b) for b in beta_s] + \
  File "/opt/anaconda3/lib/python3.7/site-packages/mpmath/functions/hypergeometric.py", line 224, in hyper
    elif q == 0: return ctx._hyp1f0(a_s[0][0], z)
  File "/opt/anaconda3/lib/python3.7/site-packages/mpmath/ctx_mp_python.py", line 1035, in f_wrapped
    retval = f(ctx, *args, **kwargs)
  File "/opt/anaconda3/lib/python3.7/site-packages/mpmath/functions/hypergeometric.py", line 271, in _hyp1f0
    return (1-z) ** (-a)
  File "<string>", line 9, in __pow__
  File "/opt/anaconda3/lib/python3.7/site-packages/mpmath/libmp/libelefun.py", line 339, in mpf_pow
    reciprocal_rnd[rnd]), -tman, prec, rnd)
  File "/opt/anaconda3/lib/python3.7/site-packages/mpmath/libmp/libmpf.py", line 1073, in mpf_pow_int
    return mpf_div(fone, inverse, prec, rnd)
  File "/opt/anaconda3/lib/python3.7/site-packages/mpmath/libmp/libmpf.py", line 960, in mpf_div
    raise ZeroDivisionError
ZeroDivisionError


我不明白这组特定值的错误原因,希望通过数值或分析积分找到解决方法。

奇怪的是,非常接近的参数值不会导致任何错误。例如,以下两组工作得很好: a1 = 1018, a2 = 1018, a3 = 917 a1 = 984, a2 = 984, a3 = 984

我要补充一点,将 a1、a2、a3 视为未知数的封闭形式解不存在。但是当用数值替换 a1、a2、a3 时,确实存在封闭形式的解决方案。问题是,根据我的需要,我需要尝试 a1、a2、a3 值的不同组合数千次。

某处存在错误。它无法尝试对此进行数值评估:

meijerg(((0, -1), ()), ((-1/2, 0), ()), 904401*exp_polar(0)/1002001)

我不确定这是否是来自符号积分例程中错误的无效表达式,或者错误是否在 mpmath 中无法对其进行评估。

该代码本身被称为缓存的假设的一部分。如果您在同一进程中重复 运行 相同的代码,则不会重新引发相同的错误。第三次 运行 我得到了答案:

In [4]: I1 = re(N(2*np.pi*a1*a2*a3*integrate(1/((a1**2 + s) * Delta(s)), (s, 0, oo))))
   ...: I2 = re(N(2*np.pi*a1*a2*a3*integrate(1/((a2**2 + s) * Delta(s)), (s, 0, oo))))
   ...: I3 = re(N(2*np.pi*a1*a2*a3*integrate(1/((a3**2 + s) * Delta(s)), (s, 0, oo))))

In [5]: I1, I2, I3
Out[5]: (4.10232880720189, 4.10232880720189, 4.3617129999554)

请注意,sympy 中对积分本身的数值评估给出了完全相同的答案,例如(使用 Integral 而不是 integrate):

In [7]: expr = 2*pi*a1*a2*a3*Integral(1/((a1**2 + s) * Delta(s)), (s, 0, oo))

In [8]: expr
Out[8]: 
             ∞                                 
             ⌠                                 
             ⎮               1                 
1905805902⋅π⋅⎮ ───────────────────────────── ds
             ⎮   ____________              2   
             ⎮ ╲╱ s + 904401 ⋅(s + 1002001)    
             ⌡                                 
             0                                 

In [9]: expr.evalf()
Out[9]: 4.10232880720189

如果您需要多次执行此操作,那么我建议使用 evalf 而不是 integrate