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: LibrariesSeverity: 3 - Normal
Priority: 5 - NormalItem Group: Inaccurate Result
Status: Patch SubmittedAssigned to: David Bateman <dbateman>
Originator Name: Pramod KumbharOriginator Email: -unavailable-
Open/Closed: OpenRelease: 3.4.3
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)

Mon 17 Mar 2014 09:15:57 PM UTC, comment #13:

Ping!

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.

D.

David Bateman <dbateman>
Project MemberIn charge of this item.
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

Marco Caliari <caliari>
Project Member
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

D.g

(file #26853)

David Bateman <dbateman>
Project MemberIn charge of this item.
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

Marco Caliari <caliari>
Project Member
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).

D.

David Bateman <dbateman>
Project MemberIn charge of this item.
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

Marco Caliari <caliari>
Project Member
Fri 26 Oct 2012 01:08:04 PM UTC, comment #5:

Couldn't this line

be simplified into

?

Marco

Marco Caliari <caliari>
Project Member
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.

D.

(file #26829)

David Bateman <dbateman>
Project MemberIn charge of this item.
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

Marco Caliari <caliari>
Project Member
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

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 MemberIn charge of this item.
Wed 24 Oct 2012 12:59:30 PM UTC, 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 Member
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))
.....
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 #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 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 4 latest changes.

    Date Changed By Updated Field Previous Value => Replaced By
    Mon 29 Oct 2012 08:48:28 PM UTCdbatemanAttached File-=>Added changeset, #26853
    Thu 25 Oct 2012 09:33:08 PM UTCdbatemanStatusNone=>Patch Submitted
    Thu 25 Oct 2012 09:20:49 PM UTCdbatemanAttached File-=>Added changeset, #26829
      Assigned toNone=>dbateman

    Back to the top


    Powered by Savane 3.1-cleanup