Implementing AbsTol for quadcc

[Rik]

9/22/17

I've been working with Nicholas on bug #52074
(https://savannah.gnu.org/bugs/?52074) to implement Matlab-compatible
functions integral2, integral3 which in Octave are wrapper functions which
call down to one of the quadrature routines quad, quadl, quadv, quadgk, or
quadcc.

The default integrator in Octave is quadcc which is a very good all-purpose
quadrature routine in that it can handle discontinuities, singularities,
NaNs, Infinite lower or upper bounds, etc.  It also happens to be written
in C++ so it is faster than quadl, quadv, and quadgk.  However, the
algorithm for quadcc only implements relative tolerance RelTol.  The
calling form is quadcc (f, a, b, reltol).  Because quadcc is not a Matlab
function, there is no requirement to have the same functional interface,
but all of the other routines use quadXX (f, a, b, abstol).  If you want to
swap back and forth between different integrators because of certain
properties of your integrand then this is inconvenient.  It also breaks
good UX conventions which favor consistency.

In addition, by only implementing RelTol there are some integrals which
become pathological.  For example, bug report #46439
(https://savannah.gnu.org/bugs/?func=detailitem&item_id=46349) shows that a
periodic function, integrated over exactly one complete period, can take an
exhaustively long time.  This code:

tic; [q] = dblquad(@(x,y) sin(2 * pi * x ) * sin(2 * pi * y),0,1,0,1, [],
'quadcc'); toc

takes 35 seconds, versus .163 seconds for quadgk.

It wasn't hard, so I made a patch to implement both absolute and relative
tolerance.  Now, the double integral above completes in .013 seconds or 12X
faster than quadgk.  I wanted to have access to both tolerances for the
algorithm so I did this

+The optional argument @var{tol} is a 1- or 2-element vector that specifies the
+desired accuracy of the result.  The first element of the vector is the
desired
+absolute tolerance, and the second element is the desired relative tolerance.
+To choose a relative test only, set the absolute tolerance to zero.  To choose
+an absolute test only, set the relative tolerance to zero.  The default
+absolute tolerance is 0 and the default relative tolerance is @math{1e^{-6}}.

This is similar to the interface for quad().

How do people feel about making this change?  Should I default it to using
[RelTol AbsTol] so existing code can run unmodified?  I'd prefer not to
since the interface would then differ from the other quad routines.  Should
there be a warning in place for two versions of Octave when the single
tolerance syntax is used informing the user that this value will now be
used for AbsTol rather than RelTol?  Should the AbsTol default to something
more sensible than 0, such as 1e-10 which is what quadgk uses?

--Rik