pythonscipynumerical-integrationquad

How to limit number of function calls in scipy.integrate.quad?


Scipy.integrate.quad() seems to make too many function calls in come cases. Here is a simple test to demonstrate:

import numpy as np
from scipy import integrate

def intgnd(x):
    p = x + x**2
    return p

x0=-1
x1=1
epsrel=0.1
epsabs=0.1
I,err,info = integrate.quad(intgnd,x0,x1,full_output=1,epsabs=epsabs,epsrel=epsrel)
print("{:.3f}, {:.3g}, {}, {}".format(I,err,info['neval'],info['last']))

The function to integrate is a second-degree polynomial, and can be integrated exactly by two-point Gaussian quadrature. quad() gets the right answer, but uses 21 points to do it. Here's the output:

0.667, 1.11e-14, 21, 1

Furthermore I only asked for an absolute error of 0.1, so I didn't even need an exact result. I have found no way to force quad() to use anything less than 21 function calls, regardless of the values of epsabs, epsrel or limit. How can I get quad() to use fewer function calls?


Solution

  • These are the x values from all 21 calls to intgnd

    In [323]: np.array(j)
    Out[323]: 
    array([ 0.        , -0.97390653,  0.97390653, -0.86506337,  0.86506337,
           -0.67940957,  0.67940957, -0.43339539,  0.43339539, -0.14887434,
            0.14887434, -0.99565716,  0.99565716, -0.93015749,  0.93015749,
           -0.78081773,  0.78081773, -0.56275713,  0.56275713, -0.29439286,
            0.29439286])
    

    With [-1,1] bounds, it appears to evaluate the function on successively narrower set of points.