如何设置 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 的幕后可能发生了很多事情，但是，评估还是在几秒钟内完成，没有任何问题。

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)

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

How to set Python scipy's nquad options to integrate near a singularity (avoid divide by zero)?

python

arguments

scipy

divide-by-zero

python-3.x