[Octave-bug-tracker] [bug #37613] Octave precision/accuracy is much lower for quadgk

[Michael Leitner]

Follow-up Comment #25, bug #37613 (project octave):

I am asking to reopen this old report instead of opening a new one, as the
relevant discussion is here:

I did some work on quadgk (bug #62468), and while doing that I wondered about
the case that sets the too-close condition. First, I do not see why it is
necessary to break off at all. If the intervals become so small that you start
to see the discretization of the numerical representation of numbers, of
course the estimate for the integral of the subinterval will become quite bad
(in a relative way). This will stay like that until the subinterval is so
small that it covers just two or even one floating-point number, at which
point the Gauss and Kronrod estimates will obviously agree, and the
subinterval will be accepted (and in my point of view, this is the correct
thing to do). More probably, however, is that the max_subintervals limit will
be reached, which will give a nice warning, and emit a sensible estimate of
the error. So I do not see why Shampine says "We break off the computation
...", and there is no explanation given. 

As it is now (seemingly, see below), if the subinterval becomes to small, its
contribution *as well as its expected contribution to the error* is set to
zero, and there is no warning. Thus, the algorithm returns as if there was no
problem, but it just cuts out those subintervals from the integration range.
This will definitely not be what Shampine meant with "We break off the
computation...", rather I would interpret this to mean to raise an error. 

*However*: in reality, the code as active since 2016 does not test this. The
line in question is 

++
if (any (abs (diff (x, [], 2) ./ max (abs (x), [], 2))) < 100 * eps1)
--

and if you look closely, you will see that the closing parenthesis of any() is
misplaced -- effectively you test whether all left and right entries of the
x-table are identical (only then the any () will not return true, that means
0, which indeed is smaller than 100 eps). So actually since 2016 this dropping
of subintervals was not applied. There have not been any reports of problems
about that, in fact I would expect to have had more reports if subintervals
where indeed dropped. 

So if nobody has an idea why it is a good idea to "break off the computation"
-- and in this case one would have to specify what this should mean, probably
not just silently dropping the subinterval -- I would just say to excise this
section of the code.


    _______________________________________________________

Reply to this item at:

  <https://savannah.gnu.org/bugs/?37613>

_______________________________________________
  Message sent via Savannah
  https://savannah.gnu.org/