python - 如何设置 Python scipy 的 nquad 选项以在奇点附近积分（避免除以零）？

Question

我正在尝试使用 Python 的 scipy 库来集成某个函数，其中涉及在c = +1. 因此，我想集成到c = 0.99，但我不知道如何设置各种选项和参数，以便集成工作。

这是一个最小的例子：

from scipy.integrate import nquad

options={'limit':5000,'epsrel':0.5e0}

results = []

for d in range(-4,5):
    f = lambda a,b,c: a**2*b**2 / (a**2 - 2*a*b*c + b**2 + d)
    temp = nquad( f, [[0,30],[0,30],[-1,0.99]], opts=[options,options,options] )
    results.append( [d, 4*10**(-8) * temp[0]] )

print(results)

我试图增加限制，但这似乎没有帮助。我也玩过epsrel价值，但无济于事。

对于它的价值，我设法在 Mathematica 中很容易地做到这一点，所以我知道这是可能的。我认为这只是我如何选择nquad选项的问题。作为参考，这是 Mathematica 的输出：

中的幕后可能发生了很多事情NIntegrate，但评估仍然在几秒钟内完成，没有任何麻烦。

标签： pythonpython-3.xscipyargumentsdivide-by-zero

Answer 1

你的函数是类型的等边双曲线

在您的情况下，奇点是发生在的垂直渐近线

所以，为了避免奇点，你可以定义你f的

f = lambda a,b,c: (a**2 * b**2) / (a**2 - 2*a*b*c + b**2 + d) \
                  if d != -(a**2 - 2*a*b*c + b**2) else 0

所以我们得到了

import numpy as np
from scipy.integrate import nquad

options={'limit':5000, 'epsrel':0.5e0}

results = []

for i, d in enumerate(np.arange(-4, 5)):
    f = lambda a,b,c: (a**2 * b**2) / (a**2 - 2*a*b*c + b**2 + d) \
                      if d != -(a**2 - 2*a*b*c + b**2) else 0
    temp = nquad( f, 
                 ranges=((0, 30), (0, 30), (-1, .99)),
                 opts=[options,options,options],
                )
    results.append( [d, 4*10**(-8) * temp[0]] )

res_arr = np.array(results)
print(res_arr)

这给出了（您可能会收到一些警告）与 Mathematica 几乎相同的结果

[[-4.          0.01405795]
 [-3.          0.01393407]
 [-2.          0.01370157]
 [-1.          0.01351541]
 [ 0.          0.01335587]
 [ 1.          0.01321535]
 [ 2.          0.01308802]
 [ 3.          0.01297009]
 [ 4.          0.01285942]]

和绘图

import matplotlib.pyplot as plt
plt.plot(res_arr[:,0], res_arr[:,1], 'k.')
plt.axvline(0, color='r', ls='--', lw=1)