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 Oct 22 13:59:36 2012  
 
Category: LibrariesSeverity: 3 - Normal
Priority: 5 - NormalItem Group: Inaccurate Result
Status: FixedAssigned to: Marco Caliari <caliari>
Originator Name: Pramod KumbharOriginator Email: -unavailable-
Open/Closed: ClosedRelease: dev
Operating System: Any

Add a New Comment (Rich MarkupRich Markup):
   

You are not logged in

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

 

(Jump to the original submission Jump to the original submission)

Tue Sep 13 04:05:40 2016, 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 Jul 11 02:43:31 2016, 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>
Project Member
Mon Jul 11 02:17:06 2016, 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>
Project Member
Thu Jun 30 17:30:37 2016, 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 Jun 30 17:21:38 2016, 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 Jun 30 15:51:21 2016, 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 ?

David

Anonymous
Thu Jun 30 14:56:09 2016, comment #18:

What about the attached patch.

(file #37633)

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

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

David

Anonymous
Thu Jun 30 13:47:09 2016, 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 Jun 30 13:07:32 2016, 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>
Project Member
Tue Dec 2 04:06:46 2014, comment #14:

Ping!

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

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

Ping!

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

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

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

Rik <rik5>
Project Administrator
Tue Oct 30 11:30:01 2012, 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.

D.

David Bateman <dbateman>
Project Member
Tue Oct 30 08:10:19 2012, 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

Marco Caliari <caliari>
Project MemberIn charge of this item.
Mon Oct 29 20:48:28 2012, 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

D.g

(file #26853)

David Bateman <dbateman>
Project Member
Mon Oct 29 08:28:05 2012, 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

Marco Caliari <caliari>
Project MemberIn charge of this item.
Fri Oct 26 20:26:55 2012, 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).

D.

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

Another thing,

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

Marco

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

Couldn't this line

be simplified into

?

Marco

Marco Caliari <caliari>
Project MemberIn charge of this item.
Thu Oct 25 21:20:49 2012, 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.

D.

(file #26829)

David Bateman <dbateman>
Project Member
Thu Oct 25 08:29:27 2012, 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

Marco Caliari <caliari>
Project MemberIn charge of this item.
Wed Oct 24 22:05:38 2012, 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

http://faculty.smu.edu/shampine/rev2vadapt.pdf

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?

D.

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

Hi,

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

Marco Caliari <caliari>
Project MemberIn charge of this item.
Mon Oct 22 13:59:36 2012, 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))
.....
break
endif

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?

Regards,
Pramod.

Anonymous

 

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

Attach File(s):
   
   
Comment:
   

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

 

Depends on the following items: None found

Items that depend on this one: None found

 

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

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

    Only project members can vote.

     

    Please enter the title of George Orwell's famous dystopian book (it's a date):

     

     

    Follow 10 latest changes.

    Date Changed By Updated Field Previous Value => Replaced By
    Tue Sep 13 04:05:40 2016rik5StatusPatch Submitted=>Fixed
      Open/ClosedOpen=>Closed
    Mon Jul 11 02:43:31 2016lachlanAttached File-=>Added bug_37613_quadgk.cset, #37795
    Mon Jul 11 02:17:06 2016lachlanAssigned todbateman=>caliari
      Release3.4.3=>dev
    Thu Jun 30 14:56:09 2016caliariAttached File-=>Added quadgk.cset, #37633
    Mon Oct 29 20:48:28 2012dbatemanAttached File-=>Added changeset, #26853
    Thu Oct 25 21:33:08 2012dbatemanStatusNone=>Patch Submitted
    Thu Oct 25 21:20:49 2012dbatemanAttached File-=>Added changeset, #26829
      Assigned toNone=>dbateman

    Back to the top


    Powered by Savane 3.1-cleanup