SymPy符号积分特定参数下报ZeroDivisionError问题求助

我使用符号积分求解以下一系列同类积分的数值结果，由于对应函数收敛速度慢，无论尝试何种相对、绝对容差设置，现有数值方法得到的结果都与符号计算结果偏差较大：

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参数组合进行求解。

内容的提问来源于stack exchange，提问作者Taiyi Wang