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

Question

I'm trying to use Python's scipy library to integrate a certain function, which involves dividing by zero when c = +1. Therefore, I want to integrate up to c = 0.99, but I cannot figure out how to set the various options and parameters, such that the integration works.

Here's a minimal example:

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)

I've tried to increase the limit, but this does not seem to help. I've also played around with the epsrel value, to no avail.

For what it's worth, I managed to do this quite easily in Mathematica, so I do know that it's possible. I assume it's just a matter of how I choose nquad's options. For reference, this is Mathematica's output:

There's probably a lot happening behind the scenes in NIntegrate, but still, the evaluation gets done in a few seconds without any trouble.

Answer 1

Your function is an equilateral hyperbola of type

the singularity in your case is the vertical asymptote that happens for

so, to avoid the singularity, you can define your f as

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

so we got

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)

that gives (you may get some warnings) almost the same results of 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]]

and plotting

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