Re: Better quadrature routine in octave

[Pedro Gonnet]

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?

> > In any case, is there any literature anywhere on how I should prepare
> > the code (formatting, comments, copyrights, etc...) before submitting it
> > to the integration package?
> >   
> Yeah the manual has a section on this.. Though the manual on the website 
> is wayyyyy out of date, so see
> 
> http://hg.savannah.gnu.org/hgweb/octave/file/tip/doc/interpreter/contrib.txi
> 
> instead.

Ok, I will have a look and patch the code up accordingly.

Cheers,
Pedro