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

[David Bateman]

Update of bug #37613 (project octave):

             Assigned to:                    None => dbateman               

    _______________________________________________________

Follow-up Comment #4:

If there is an infinite interval, then h0 is always infinite and the test 


if (any (abs (diff (trans (subs), [], 2) / h0) < 100 * myeps))


is always going to suceed and the quadrature exit immediately. I suppose
you're right Marco in that matlab must test at the abscissa in the
untransformed subintervals rather than do a single overall test on the bounds
of the subintervals. Though Shampine doesn't make it clear why the test *MUST*
be done in the untransformed intervals I supposed its safer to do so.

I don't think the changes necessary will be that major. We just need to move
the test into __quadgk_eval__ where the abscissa are available. As we would
then be dealing with a matrix of abscissa rather than a vector of sub interval
end points the test would either need to use bsxfun or the new broadcasting
code.

I propose the attached changeset, though a little testing before committing it
might be good idea. In any case with this changeset I now see


octave:24> f = @(x)x.^5.*exp(-x).*sin(x); 
octave:25> [q,errbnd] = quadgk(f,0,inf,'RelTol',1e-8,'AbsTol',1e-12) 
q = -15.000
errbnd =  1.0290e-08


So it appears to work correctly.

D.



(file #26829)
    _______________________________________________________

Additional Item Attachment:

File name: changeset                      Size:3 KB


    _______________________________________________________

Reply to this item at:

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

_______________________________________________
  Message sent via/by Savannah
  http://savannah.gnu.org/