BUG: integrate: ensure evaluates the function at least once; also, ad…

…d a relative tolerance in stopping condition

Also, use warnings instead of printing to stdout in quadrature

scipy/integrate/quadrature.py

@@ -8,6 +8,10 @@
      asarray, real, trapz, arange, empty
 import numpy as np
 import math
+import warnings
+
+class AccuracyWarning(Warning):
+    pass
 
 def fixed_quad(func,a,b,args=(),n=5):
     """
@@ -99,7 +103,8 @@ def vfunc(x):
             return output
     return vfunc
 
-def quadrature(func,a,b,args=(),tol=1.49e-8,maxiter=50, vec_func=True):
+def quadrature(func, a, b, args=(), tol=1.49e-8, rtol=1.49e-8, maxiter=50,
+               vec_func=True):
     """
     Compute a definite integral using fixed-tolerance Gaussian quadrature.
 
@@ -116,9 +121,9 @@ def quadrature(func,a,b,args=(),tol=1.49e-8,maxiter=50, vec_func=True):
         Upper limit of integration.
     args : tuple, optional
         Extra arguments to pass to function.
-    tol : float, optional
+    tol, rol : float, optional
         Iteration stops when error between last two iterates is less than
-        tolerance.
+        `tol` OR the relative change is less than `rtol`.
     maxiter : int, optional
         Maximum number of iterations.
     vec_func : bool, optional
@@ -147,17 +152,20 @@ def quadrature(func,a,b,args=(),tol=1.49e-8,maxiter=50, vec_func=True):
     odeint: ODE integrator
 
     """
-    err = 100.0
-    val = err
-    n = 1
     vfunc = vectorize1(func, args, vec_func=vec_func)
-    while (err > tol) and (n < maxiter):
+    val = np.inf
+    err = np.inf
+    for n in xrange(1, maxiter+1):
         newval = fixed_quad(vfunc, a, b, (), n)[0]
         err = abs(newval-val)
         val = newval
-        n = n + 1
-    if n == maxiter:
-        print "maxiter (%d) exceeded. Latest difference = %e" % (n,err)
+
+        if err < tol or err < rtol*abs(val):
+            break
+    else:
+        warnings.warn(
+            "maxiter (%d) exceeded. Latest difference = %e" % (maxiter, err),
+            AccuracyWarning)
     return val, err
 
 def tupleset(t, i, value):

scipy/integrate/tests/test_quadrature.py

@@ -1,8 +1,7 @@
-
 import numpy
 from numpy import cos, sin, pi
 from numpy.testing import TestCase, run_module_suite, assert_equal, \
-    assert_almost_equal
+    assert_almost_equal, assert_allclose
 
 from scipy.integrate import quadrature, romberg, romb, newton_cotes
 
@@ -18,6 +17,13 @@ def myfunc(x,n,z):       # Bessel function integrand
         table_val = 0.30614353532540296487
         assert_almost_equal(val, table_val, decimal=7)
 
+    def test_quadrature_rtol(self):
+        def myfunc(x,n,z):       # Bessel function integrand
+            return 1e90 * cos(n*x-z*sin(x))/pi
+        val, err = quadrature(myfunc,0,pi,(2,1.8),rtol=1e-10)
+        table_val = 1e90 * 0.30614353532540296487
+        assert_allclose(val, table_val, rtol=1e-10)
+
     def test_romberg(self):
         # Typical function with two extra arguments:
         def myfunc(x, n, z):       # Bessel function integrand