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

如何解决如何设置 Python scipy 的 nquad 选项以在奇点附近进行积分避免除以零？

我正在尝试使用 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,[-1,0.99]],opts=[options,options,options] )
    results.append( [d,4*10**(-8) * temp[0]] )

print(results)

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

不管怎样，我在 Mathematica 中很容易地做到了这一点，所以我知道这是可能的。我认为这只是我如何选择 nquad 的选项的问题。作为参考，这是 Mathematica 的输出：

NIntegrate 的幕后可能发生了很多事情，但仍然可以在几秒钟内轻松完成评估。

解决方法

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

你的情况的奇点是发生在

的垂直渐近线

因此，为了避免奇点，您可以将 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,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,(-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)