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!

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)

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"?
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)?

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.

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.

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

What about the attached patch.
(file #37633)

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

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

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

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

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

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

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.

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

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)

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

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 subinterval 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.

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

Couldn't this line
be simplified into
?
Marco

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)

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

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.

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

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',1e8,'AbsTol',1e12)
q = 1.499999999999836e+01
errbnd = 9.438576028474561e09
But in the Octave, I get
debug> [q,errbnd] = quadgk(f,0,inf,'RelTol',1e8,'AbsTol',1e12)
warning: quadgk: Error tolerance not met. Estimated error 0.0329476
q = 1.49999941897813e+01
errbnd = 3.29475856720323e02
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/octave3.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.
