# 2D integrals in Python with specified points of interest

Note: This is my first question on stackexchange; please tell me if I'm doing something incorrectly.

I am trying to calculate a series of a 2D integrals in Python with an integrand that has several highly localized peaks. I've been using scipy.integrate.dblquad. However, repeatedly evaluating this integral for a series of integrands with successively more localized peaks is taking a long, long time.

I'm hoping I can speed things up by using a 2D adaptive integration method which allows me to specify an approximate location for these peaks, similar to the points optional input in scipy.integrate.quad. Does anyone know of an easy way to accomplish this?

Denoting my integrand as some scalar function f(x,y), what I currently have looks something like:

import numpy as np
from scipy import integrate
from joblib import Parallel, delayed
import multiprocessing

num_cores = multiprocessing.cpu_count()

#Values V for which I evaluate f(x,y,V)
V_list = np.arrange(0,5,0.1)

#Setting up a function to evaluate the integrals at value V
def integration_process(V):
integrand = lambda x, y: f[x,y,V]
return output

#Using all available cores to compute integrals in parallel
results = Parallel(n_jobs=num_cores,verbose=10)(delayed(integration_process)(V) for V in V_list)


I could probably get the points of interest with scipy.optimize.root_scalar operating on part of the integrand, but I can't find an adaptive integration method in scipy that will let me use that information to help determine the initial integration points.

Digging through the source code of the scipy integration routines I found most of them (all i looked at) are using QUADPACK under the hood, which provides the following functionality:

nquad calls routines from the FORTRAN library QUADPACK. This section provides details on the conditions for each routine to be called and a short description of each routine. The routine called depends on weight, points and the integration limits a and b.

================  ==============  ==========  =====================
QUADPACK routine  weight          points      infinite bounds
================  ==============  ==========  =====================
qagse             None            No          No
qagie             None            No          Yes
qagpe             None            Yes         No
qawoe             'sin', 'cos'    No          No
qawfe             'sin', 'cos'    No          either a or b
qawse             'alg*'          No          No
qawce             'cauchy'        No          No
================  ==============  ==========  =====================


I see that scipy.integrate.dblquad is just a wrapper for nquad, but doesn't offer any way of passing along extra options:

    def ranges0(*args):
return [qfun(args[1], args[0]) if callable(qfun) else qfun,
rfun(args[1], args[0]) if callable(rfun) else rfun]

def ranges1(*args):
return [gfun(args[0]) if callable(gfun) else gfun,
hfun(args[0]) if callable(hfun) else hfun]

ranges = [ranges0, ranges1, [a, b]]
opts={"epsabs": epsabs, "epsrel": epsrel})


but nquad does mention supporting "points" via opts

opts: iterable object or dict, optional

Options to be passed to quad. May be empty, a dict, or a sequence of dicts or functions that return a dict. If empty, the default options from scipy.integrate.quad are used. If a dict, the same options are used for all levels of integraion. If a sequence, then each element of the sequence corresponds to a particular integration. e.g., opts[0] corresponds to integration over x0, and so on. If a callable, the signature must be the same as for ranges. The available options together with their default values are:

• epsabs = 1.49e-08
• epsrel = 1.49e-08
• limit = 50
• points = None
• weight = None
• wvar = None
• wopts = None

For more information on these options, see quad.

so it should be possible to use that directly:

integrate.nquad(integrand, [[0.001, 10], [-10, 10]], opts=[{'points': x_points}, {'points': y_points}])


You could even use different epsabs or epsrel, limits, weights for the x and y axis (if that makes sense for your particular function).

## Limitations

The way QUADPACK works (it's only a 1-D integration engine) inside scipy.integrate, I believe this effectively treats the entire lines at the given x and y points special. I.e. if you have a singularity at the point (x,y) = (1,1), this effectively tells QUADPACK to treat the lines x=1 as special, and y=1 as special. That should still be fine I think. I also have no idea if this actually speeds anything up for near singular cases.

From my understanding of the scipy code, there wouldn't be any way to handle lines of (near) singularity along an edge that isn't parallel with the x or y axis.

At the core with anything adaptive, I would expect it to become a lot more expensive as you make your function more evil. You could adapt the epsabs and/or epsrel depending on what you find acceptable for your application. Alternatively, you could use fixed integration grid of some sort and just accept the loss of precision if you need an estimate within a fixed time.