[Octave-bug-tracker] [bug #59403] integral2 and quad2d functions do not accurately document form of integrand function F

[Nicholas Jankowski]

Follow-up Comment #11, bug#59403 (group octave):

I didn't really give this much attention a few years ago.  (trying to
graduate, etc.)

reading through things to see where a doc change might be made, I'm think that
quad2d may actually be doing something wrong when it's passed a
not-properly-vectorized function (as opposed to it just doing a proper
integral on the non-elementwise function it's given). for this f, only quad2d
and the integral2 and 3 functions calling quad2d are generating different
numerical outputs without erroring. generally, it's important that integration
query points be independent (otherwise the summation loses meaning), and
sinc(x)sinc(y) still does that when y is a scalar. integral2 and 3 select
between iterated (using quadcc), and quad2d. Iterated produces a correct
result.  quad2d doesn't. 

The comments below basically capture the situation. When I made the 1st stab
at integral2, I retained the input requirements from dblquad (vector x, scalar
y) as it was really just a wrapper.  then we (well david) improved it to
actually do a decent job, creating quad2d, too, but we mistakenly kept those
same input requirements for f in the doc.  As rik saw, stepping through
quad2d, at the end of the tensorproduct subfunction it calls f(x,y) after
expanding both x and y to  2D matrices like they're doing a manual 'meshgrid'.
That appears in fact to be a feature of the Shampine algorithm our quad2d (and
matlab's according to their help doc) is based on. This sort of requires
elementwise operation so it can do vectorized expansion for speed. if the
elements aren't independent, it loses meaning. iterated methods sometimes get
away with odd functions, but by permitting non-elementwise operations, we're
really violating the intent of the Shampine method. most of the integrators
have options to turn off vectorized solving specifically to navigate around
this problem.

right now 1d, 2d, and 3d integrators wind up failing with obscure errors for
malformed f's, catching them would be better.  And apparently quad2d may have
been silently producing erroneous results since octave 4. from jwe's comment
#2, I'm guessing we don't have the ability to throw an upfront error? that
would be ideal. the only thing I could think of is something crude like
assert(vectorize(f),f).  but even vectorize says not to use vectorize. and
that won't catch much anyway.

If the best we can do is a doc change, i think it needs to state tighter input
requirements, specifically as elementwise.  "inputs need to have x form or
numerical results can be unpredictable. if you can't then use
'vectorized=false'" that way at least people are warned, with a workaround,
and silent failures like comment #0 can be waved away with 'we told you so'. 
I'll generate some text and see how it reads.


    _______________________________________________________

Reply to this item at:

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

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