[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
Re: Better quadrature routine in octave
From: |
Pedro Gonnet |
Subject: |
Re: Better quadrature routine in octave |
Date: |
Tue, 30 Mar 2010 23:03:38 +0100 |
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
- Better quadrature routine in octave, Pedro Gonnet, 2010/03/30
- Re: Better quadrature routine in octave, David Bateman, 2010/03/30
- Message not available
- Re: Better quadrature routine in octave, Pedro Gonnet, 2010/03/31
- Re: Better quadrature routine in octave, Jaroslav Hajek, 2010/03/31
- Re: Better quadrature routine in octave, Michael D Godfrey, 2010/03/31
- Re: Better quadrature routine in octave, Pedro Gonnet, 2010/03/31
- Re: Better quadrature routine in octave, Michael D Godfrey, 2010/03/31
- Re: Better quadrature routine in octave, Pedro Gonnet, 2010/03/31