bugGNU Octave - Bugs: bug #37613, Octave precision/accuracy is much...


bug #37613: Octave precision/accuracy is much lower for quadgk

Submitted by:  None
Submitted on:  Mon 22 Oct 2012 01:59:36 PM UTC  
Category:  Libraries Severity:  3 - Normal
Priority:  5 - Normal Item Group:  Inaccurate Result
Status:  Fixed Assigned to:  Marco Caliari <caliari>
Originator Name:  Pramod Kumbhar Originator Email:  -email is unavailable-
Open/Closed:  Closed Release:  dev
Operating System:  Any

Add a New Comment(Rich Markup)

You are not logged in

Please log in, so followups can be emailed to you.


( Jump to the original submission)

Tue 13 Sep 2016 04:05:40 AM UTC, comment #24:

I finally pushed Marco's patch by way of Lachlan here (http://hg.savannah.gnu.org/hgweb/octave/rev/ed1722f70fad). Closing the bug from 2012!

Rik <rik5>
Project Administrator
Mon 11 Jul 2016 02:43:31 AM UTC, comment #23:

I forgot to provide an updated patch with the last comment. Here it is. (It keeps Marco as the author.)

Marco's is fine without my comments, so take your pick.

(file #37795)

Lachlan Andrew <lachlan>
Mon 11 Jul 2016 02:17:06 AM UTC, comment #22:

Some trivial feedback on the patch:

  • The new assert in the BIST is much longer than 80 characters. Could it rewritten to fit under 79 characters? Perhaps defining

would help.

  • Is it worth retaining the comment "Shampine suggests 100*eps"?
  • The code

hides the fact that this is an early return. How about

  • Should "tooclose" be "too_close"?
  • I don't like shadowing the function eps by a variable. I know that was in the existing code, but it is also done in the new _quadgk_eval_. How about calling it eps1, since it is eps (1)?
Lachlan Andrew <lachlan>
Thu 30 Jun 2016 05:30:37 PM UTC, comment #21:

Less than 1% change basically means no difference. That small a time difference could be caused by various background processes waking up on your computer, etc.

Seems like it is worth committing the patch.

Rik <rik5>
Project Administrator
Thu 30 Jun 2016 05:21:38 PM UTC, comment #20:

Actually, I got an average of 0.45156 for the 24 tests of the new version vs. 0.45522 for the 23 tests of the old version.

Marco Caliari <caliari>
Project MemberIn charge of this item.
Thu 30 Jun 2016 03:51:21 PM UTC, comment #19:

I'm too far out of Octave developpement to be able to give you a definitive yes or no. However, my original objective to this proposed change what the requiring the call to @trans for all of the tests will make it slower. It would be interesting to know how much slower. What does

tic; test guadgk; toc

give before and after your change for exemple ?


Thu 30 Jun 2016 02:56:09 PM UTC, comment #18:

What about the attached patch.

(file #37633)

Marco Caliari <caliari>
Project MemberIn charge of this item.
Thu 30 Jun 2016 02:06:36 PM UTC, comment #17:

In any case Marco's comments 8 and 10 identified problems in my patch that should be fixed before it is applied


Thu 30 Jun 2016 01:47:09 PM UTC, comment #16:

The issue is not resolved and David's path does not apply anymore. I can try to submit a new patch.

Marco Caliari <caliari>
Project MemberIn charge of this item.
Thu 30 Jun 2016 01:07:32 PM UTC, comment #15:

Ping, ping, ping!

David, would you like this patch to get in to 4.2.0, or has this issue been resolved?

Lachlan Andrew <lachlan>
Tue 02 Dec 2014 04:06:46 AM UTC, comment #14:


David, do you want to commit the patch you prepared for this bug report?

Rik <rik5>
Project Administrator
Mon 17 Mar 2014 09:15:57 PM UTC, comment #13:


David, do you want to commit the patch you prepared for this bug report?

Rik <rik5>
Project Administrator
Mon 21 Oct 2013 04:38:23 AM UTC, comment #12:

David, do you want to go ahead and apply your changeset from last year?

Rik <rik5>
Project Administrator
Tue 30 Oct 2012 11:30:01 AM UTC, comment #11:

I don't think it is "equal to the old one". Please note I changed the variable "x" in _quadgk_eval_ to be "t" and made the new variable "x" be the abscissae in the untransformed space.

Yes it should be the maximum of the absolute value of x.


David Bateman <dbateman>
Project Member
Tue 30 Oct 2012 08:10:19 AM UTC, comment #10:

I agree with your comment, but your new patch is equal to the old one.
By the way, shouldn't you take the maximum of the absolute value of x?


Marco Caliari <caliari>
Project MemberIn charge of this item.
Mon 29 Oct 2012 08:48:28 PM UTC, comment #9:

Opps, yes it should be 'transformed' to the 'untransformed' form with @trans. However, thinking about it, that going to be fairly expensive if we test the distance between all of the abscissae. Also it seems difficult to understand what Shampine meant by

and my previous test was wrong for yet another reason; It trying to give a "relative value" I effectively always renormalized the distance between abscissae by the total length between the abscissae (ie 2) so that the value returned in the test would be 0.021 (ie the distance between the 1st and 2nd abscissae divided by 2) so comparing to 100 * eps would always return false.

Thinking about it Shampine probably meant to compare "DeltaX / X < 100 * eps". As the test is relatively artifical we can get the same effective result by comparing between any two abscissae in the interval and so a test like

Should do the right thing. Here is a modified changset


(file #26853)

David Bateman <dbateman>
Project Member
Mon 29 Oct 2012 08:28:05 AM UTC, comment #8:

If 'untransformed' is (0,Inf) and 'transformed' is (0,1), according to Shampine the check should be done in (0,Inf), i.e. with trans(abscissa). So, we should pass @trans to @__quadgk_eval__


Marco Caliari <caliari>
Project MemberIn charge of this item.
Fri 26 Oct 2012 08:26:55 PM UTC, comment #7:

Yes I think you're right the simplification will give the same results much faster.

In fact I'm no longer testing the minimum relative size of an interval, but rather the minimum distance between abscissa, since as you said the "arguments" mentioned by Shampine probably mean the distance between abscissa or evaluation points rather than the sub-interval length.

The test is also after the conversion from the transformed interval to the real evaluation points so I'm testing in the interval (a, b).


David Bateman <dbateman>
Project Member
Fri 26 Oct 2012 03:25:49 PM UTC, comment #6:

Another thing,

are you not testing now in the trasformed interval, i.e. [0,1]?


Marco Caliari <caliari>
Project MemberIn charge of this item.
Fri 26 Oct 2012 01:08:04 PM UTC, comment #5:

Couldn't this line

be simplified into



Marco Caliari <caliari>
Project MemberIn charge of this item.
Thu 25 Oct 2012 09:20:49 PM UTC, comment #4:

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

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

So it appears to work correctly.


(file #26829)

David Bateman <dbateman>
Project Member
Thu 25 Oct 2012 08:29:27 AM UTC, comment #3:

Dear David,

maybe the key is in the word "arguments": you took the length of subintervals as arguments (and in the case [0,inf) there is always a subinterval of infinite length), but in the paper the arguments are the evaluation points for the integrand (see, for instance, the third line in p.4). Since inf is never an evaluation point, the distance between "arguments" is never infinite.
This would require some work on the code...


Marco Caliari <caliari>
Project MemberIn charge of this item.
Wed 24 Oct 2012 10:05:38 PM UTC, comment #2:

In fact if you replaced h0 with h you'd also have to replace trans(subs) with subs. As for the sense of the comment, I believe I based it on the statement in section section of the Shampine paper

The paper itself can be found at


and I didn't have any other reasons for this comment. However changing line 322 of quadgk to read

then running "test quadgk" passes all of the tests. These tests include tests for infinite bounds and singular integrands and so probably give a pretty clear indication that such a change would be acceptable, even if Shampine himself didn't recommend it. Perhaps Mathworks modified this behavior after the inclusion of quadgk in Matlab?


David Bateman <dbateman>
Project Member
Wed 24 Oct 2012 12:59:30 PM UTC, comment #1:


instead of removing the check, I think that h0 (length of untransformed interval) should be replaced by h (length of transformed interval). But the comment relative to this check explicitly speaks of "untransformed subintervals". Maybe David had some counterexample in mind.


Marco Caliari <caliari>
Project MemberIn charge of this item.
Mon 22 Oct 2012 01:59:36 PM UTC, original submission:

Hello All,

I am new to Octave and running Matlab code on Octave 3.4.3. I am getting incorrect results for my simulations, basically because of lower precision/accuracy of Octave quadrature routine quadgk.

For example, say we have

f = @(x)x.^5.*exp(-x).*sin(x);

In Matlab I get:

K>> [q,errbnd] = quadgk(f,0,inf,'RelTol',1e-8,'AbsTol',1e-12)

q = -1.499999999999836e+01
errbnd = 9.438576028474561e-09

But in the Octave, I get

debug> [q,errbnd] = quadgk(f,0,inf,'RelTol',1e-8,'AbsTol',1e-12)
warning: quadgk: Error tolerance not met. Estimated error 0.0329476
q = -1.49999941897813e+01
errbnd = 3.29475856720323e-02

From above simple test, I see Octave accuracy is very low.

To find the cause of this, I loooked into the implementation of quadgk routine in Octave (here:http://fossies.org/dox/octave-3.6.3/quadgk_8m_source.html).

The problem is, when one of the limit of integral is Inf, current implementation don't divide subintervals further due to below condition (line no 322):

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

When a or b is Inf, h0 is Inf and hence code returns results with very low accuracy.

I am not a mathematician and not sure if this is expected behaviour. But simply removing above check we get more accurate results similar to Matlab.

Can someone please help in this issue? Suggestions for getting more accurate results?




(Note: upload size limit is set to 16384 kB, after insertion of the required escape characters.)

Attach Files:

Attached Files
file #37795:  bug_37613_quadgk.cset added by lachlan (4KiB - text/x-diff)
file #37633:  quadgk.cset added by caliari (3KiB - application/x-tex)
file #26853:  changeset added by dbateman (3KiB - application/octet-stream)
file #26829:  changeset added by dbateman (3KiB - application/octet-stream)


Depends on the following items: None found

Items that depend on this one: None found


Carbon-Copy List
  • -email is unavailable- added by lachlan (Posted a comment)
  • -email is unavailable- added by rik5 (Posted a comment)
  • -email is unavailable- added by dbateman (Posted a comment)
  • -email is unavailable- added by caliari (Posted a comment)
  • -email is unavailable- added by None (Submitted the item)

    Do you think this task is very important?
    If so, you can add your encouragement to it.
    This task has 0 encouragements so far.

    Only project members can vote.




    Follow 10 latest changes.

    Date Changed by Updated Field Previous Value => Replaced by
    2016-09-13 rik5 StatusPatch Submitted => Fixed
        Open/ClosedOpen => Closed
    2016-07-11 lachlan Attached File- => Added bug_37613_quadgk.cset, #37795
    2016-07-11 lachlan Assigned todbateman => caliari
        Release3.4.3 => dev
    2016-06-30 caliari Attached File- => Added quadgk.cset, #37633
    2012-10-29 dbateman Attached File- => Added changeset, #26853
    2012-10-25 dbateman StatusNone => Patch Submitted
    2012-10-25 dbateman Attached File- => Added changeset, #26829
        Assigned toNone => dbateman

    Back to the top

    Powered by Savane 3.3