[SciPy-Dev] severe bug in scipy.integrate.quadrature?

[Thomas Robitaille]

Hi,

I think there is a mistake in scipy.integrate.quadrature that can result in false convergence of the integral. The code in question is the following:

    err = 100.
    val = err
    n = 1
    vfunc = vectorize1(func, args, vec_func=vec_func)
    while (err > tol) and (n < maxiter):
        newval = fixed_quad(vfunc, a, b, (), n)[0]
        err = abs(newval-val)
        val = newval
        n = n + 1

The error inside the loop is the absolute error, so it looks like the tolerance to specify is the absolute tolerance, not relative. However, the original err is set to 100, so if I compute an integral that has a result of 1e20, I might want a tolerance of say 1e15, in which case the loop will be skipped, because err is not greater than tol. In its present form, I think that err should be initially set to infinity.

In addition, I think that the tolerance to specify should be the relative one, not the absolute, because most of the time, I know I want my integrals calculated to a precision of say 1e-5, but I don't know in advance the answer! I think a relative tolerance is what users expect at any rate. 

Can any experts comment on this?

Thanks,

Thomas