1

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?

Marvin Germain
  • 13
  • 2
  • Since you can read the docs as well as I can, I probably shouldn't bother to comment. Do you understand the full output? Have you looked at other integration methods? I see gausian in a couple of descriptions. – hpaulj May 18 '21 at 00:51
  • From the docs, 0.667 is the result of the integration, which I checked by hand and is correct to 3digits. The 1.11e-14 is the estimated absolute error, which is close enough to zero for me. "nval" in the number of function evaluations, and "last" is the number of subintervals, which does not go any lower than 1. – Marvin Germain May 18 '21 at 15:58
  • I agree that a possible solution is just don't use quad. But the point of my question was in case quad is fine, and I'm doing something wrong. – Marvin Germain May 18 '21 at 16:08
  • Why should there be anything wrong with what you are doing. The task is passed some "optimized" Fortran/C package. We don't know the details of how that works (I don't know anything, for that matter). Short of studying the code, the best we can do is experiment. Maybe there's isn't anything we can do to reduce the number of calls. – hpaulj May 18 '21 at 16:15
  • It's not hard to write your own crude numeric integration. We all learned the basic trapezoidal rule. "simpsons" is another name that comes to mind, though it's been years since I studied this stuff. – hpaulj May 18 '21 at 16:19
  • I have not tried to study the code, only the docs which are terse. My experiments say the minimum number of calls is 21. But I didn't know if there was a trick or if someone had faced this problem before. Naturally I want to check that before abandoning a mature package and writing my own. So just now I implemented a Gaussian 5-point quadrature. It runs 4x faster and the numbers (in my real project, not the toy example) agree to 5 decimal places. – Marvin Germain May 18 '21 at 16:52

1 Answers1

1

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.

hpaulj
  • 221,503
  • 14
  • 230
  • 353
  • I think your comment that there may not be any way to reduce the number of function calls is the actual answer. Thank you for taking the time to think about the problem. – Marvin Germain May 18 '21 at 16:55