Error using SymPy symbolic integration for a specific set of parameters

I used symbolic integration to find the numerical values of the following series of (similar) integrals, because the functions appear to be slow convergent, and numerical methods so far result in values quite different from the symbolic calculation (no matter what relative and absolute tolerance I tried):

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

The integration works for most values of a1, a2, a3. However, by mere chance, I found a set of values that results in errors : a1 = 1001, a2 = 1001, a3 = 951. With the following error:

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


I do not understand the cause for error with this specific set of values, and would like to find a way to work around it, either through numerical or analytical integration.

Curiously, parameter values that are quite close does not cause any errors. For example, the following two sets work perfectly fine: a1 = 1018, a2 = 1018, a3 = 917 a1 = 984, a2 = 984, a3 = 984

I would add that closed form solution treating a1, a2, a3 as unknowns do not exist. But closed form solutions do exist when replacing a1, a2, a3 with numerical values. The problem being that for my need I need to try different combinations of a1, a2, a3 values thousands of times.

Answer

There is a bug somewhere. It fails trying to numerically evaluate this:

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

I’m not sure if that’s an invalid expression coming from a bug in the symbolic integration routine or if the bug is in mpmath for not being able to evaluate it.

That code is itself called as part of the assumptions which are cached. The same error does not get reraised if you run the same code repeatedly within the same process. On the third run I get the answers:

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)

Note that numerical evaluation of the integral itself in sympy gives the exact same answers e.g. (use Integral rather than 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

If you need to do this many times then I would suggest using evalf rather than integrate.