Ошибка при использовании символической интеграции SymPy для определенного набора параметров. ⇐ Python

[Гость]

Я использовал символьное интегрирование, чтобы найти числовые значения следующих серий (похожих) интегралов, поскольку функции кажутся медленно сходящимися, а численные методы до сих пор приводят к значениям, сильно отличающимся от символьных вычислений (независимо от того, какие относительные и абсолютную толерантность я пробовала):

Код: Выделить всё

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 
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 "", 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.


Источник: https://stackoverflow.com/questions/690 ... parameters