Re: Better quadrature routine in octave

[Jaroslav Hajek]

On Wed, Mar 31, 2010 at 12:03 AM, Pedro Gonnet <address@hidden> wrote:
>
> Hi David,
>
>> If the test of correct or incorrect is just whether the expected
>> absolute and relative tolerance is met the fact that quadgk is not that
>> great isn't a surprise as quadgk uses a 15 point guass rule compared to
>> a 7 point rule to identify sub-intervals that it considered are
>> converged.. It seems your code is a bit more intelligent... In most
>> cases I've worked with (as an engineer rather than a mathematician) this
>> distinct between correct an incorrect is rarely important as the
>> convergence is still "good enough".. Octave uses the Quadpack function
>> xQAPI and xQAGP rather than DQAGS to allow treatment of user supplied
>> signularities.
>
> Actually, the degree of the underlying quadrature rule doesn't really
> matter that much. If you've got a good error estimate, then even
> Simpson's rule will integrate anything.
>
> Regarding the "good enough", in the tests I used a strict pass-fail
> criteria, i.e. I expected to get the number of digits I asked for. I
> should note, however, that when an algorithm failed it usually did so
> quite spectacularly. There aren't any details in the TOMS paper, but I
> have a review paper submitted to ACM Computing Surveys
> (http://arxiv.org/abs/1003.4629) which has some more detailed results.
>
> This is, however, academic, and the question of what reliability should
> be is a more philosophical issue on which probably everybody has their
> own opinion and over which I gladly defer to whatever the consensus is
> amongst Octave developers :)
>
>> It is true that in the core of Octave accuracy is preferred over speed,
>> and yes I consider that replacing quad would be a good idea, mainly
>> because the quadpack code being in F77 is not recursive and we can't
>> nest the quadratures. I suppose your code might be a good code to use
>> for a replacement, though I don't think I'm the one to convince of this,
>> but you'd have to wrap the code in a bit of code that sub-divides at the
>> singularities before calling your cqaud code on each sub-interval.
>
> That shouldn't be too difficult: The code works with an explicit heap of
> intervals so I'd just have to initialize the heap not with one, but with
> all intervals.
>
>> Writing it in C++ using the oct-file interface might be a good thing as
>> well.
>
> This should not be a problem as I was planning on writing a C-language
> version for the GNU Scientific Library anyway. I guess there is a
> documentation somewhere for using this interface?
>

I advise you not to do it until it becomes apparent that the m-file
version can't be optimized to provide satisfactory performance.
Having algorithms in m-files has certain advantages, in particular
that it is easier for people to play with them.

If you come up with a version that can replace quad (if I understood
correctly, it suffices to add handling of singular points)
and adapt it somewhat to the coding standards (at least the major
points, i.e. consistent comments, specific block ends), I can help you
with optimizations, should any be needed.

regards

-- 
RNDr. Jaroslav Hajek, PhD
computing expert & GNU Octave developer
Aeronautical Research and Test Institute (VZLU)
Prague, Czech Republic
url: www.highegg.matfyz.cz