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