Re: Better quadrature routine in octave

[David Bateman]

Please keep responses on the list

Pedro Gonnet wrote:Hi David,
> I've been asked for this comparison before, so I've already got the
> results (attached). It's essentially the same data as in the paper, but
> with quadgk too. Since quadgk uses a similar error estimate as DQAGS,
> the results of both algorithms are in fact quite similar.
>   
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.

> I guess I could add the code "as-is" to the integration package, but I
> was wondering if you were interested in replacing the default
> integrator, i.e. whatever is called with "quad". Since most users call
> this routine without thinking, my guess is that reliability should trump
> efficiency. Users who know a bit more about their integrand could then
> go to potentially less reliable routines if they need more speed.
>   
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. 
Writing it in C++ using the oct-file interface might be a good thing as 
well.

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

> Cheers and thanks,
> Pedro
>
>
>   


--
David Bateman                                address@hidden
35 rue Gambetta                              +33 1 46 04 02 18 (Home)
92100 Boulogne-Billancourt FRANCE            +33 6 72 01 06 33 (Mob)