对一组特定参数使用 SymPy 符号集成时出错

Question

我使用符号积分来求以下一系列（相似）积分的数值，因为这些函数似乎收敛缓慢，并且迄今为止的数值方法导致的值与符号计算完全不同（无论相对和我尝试过的绝对容忍度）：

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 值的不同组合数千次。

Answer 1

某处存在错误。尝试对此进行数值评估失败：

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。