对一组特定参数使用 SymPy 符号积分时出错
Error using SymPy symbolic integration for a specific set of parameters
我使用符号积分求出以下一系列（类似的）积分的数值，因为这些函数收敛速度较慢，而且迄今为止的数值方法得出的值与符号计算的值大不相同（无论我试过什么相对和绝对容忍度）：
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 值的不同组合数千次。
某处存在错误。它无法尝试对此进行数值评估：
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。
我使用符号积分求出以下一系列（类似的）积分的数值，因为这些函数收敛速度较慢，而且迄今为止的数值方法得出的值与符号计算的值大不相同（无论我试过什么相对和绝对容忍度）：
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 值的不同组合数千次。
某处存在错误。它无法尝试对此进行数值评估：
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。