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?

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?

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.

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

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)

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

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

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

Fri 26 Oct 2012 01:08:04 PM UTC, comment #5:
Couldn't this line
be simplified into
?
Marco

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)

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

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.

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

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',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.
