[rfc] Quadrature library for numerical integration of definite integrals
Hi, I have put together a quadrature library, for the numerical evaluation of definite integrals. I put this together to meet some requirements that I had. I would like to get comments on whether the scope, algorithms and approach used are of interest within boost. The code is based on QUADPACK's QAG and QAG global adaptive routines, using Kronrod-Gauss or Recursive Monotone Stable integration "kernels". This is not just a translation of the QUADPACK code, but tries to take a (reasonably) modern C++ approach. The RMS kernel is based on TOMS 691. Error estimation, choice of integration kernel, algorithm performance statistics, series limit acceleration algorithm, etc are all plugable. For the less brave, I would appreciate any comments you might have on a quadrature library in general. For the braver, the code can be found in the boost vault at http://www.boost-consulting.com/vault/index.php?direction=0&order=&directory=Math%20-%20Numerics& The jamfile expects BOOST_ROOT to be defined. The code use the math toolkit from the same place in the vault. The location of the toolkit can be specified using the BOOST_MATH_TOOLKIT_ROOT environment variable, and by default is looked for in a directory boost-math-toolkit, parallel to BOOST_ROOT (ie in BOOST_ROOT/..) The code has been tested against QUADPACK and includes (very) rudimentary documentation, with examples. The code has been compiled only on windows, using MSVC 8 and gcc 3.4.4. Hugo
Minor point, please package tar files under a top-level directory - that's what untar-ers expect.
-----Original Message----- From: boost-bounces@lists.boost.org [mailto:boost-bounces@lists.boost.org] On Behalf Of Hugo Duncan Sent: 05 December 2007 23:01 To: boost@lists.boost.org Subject: [boost] [rfc] Quadrature library for numerical integration ofdefinite integrals
Hi,
I have put together a quadrature library, for the numerical evaluation of definite integrals.
I put this together to meet some requirements that I had. I would like to get comments on whether the scope, algorithms and approach used are of interest within boost.
The code is based on QUADPACK's QAG and QAG global adaptive routines, using Kronrod-Gauss or Recursive Monotone Stable integration "kernels".
This is not just a translation of the QUADPACK code, but tries to take a (reasonably) modern C++ approach. The RMS kernel is based on TOMS 691.
Error estimation, choice of integration kernel, algorithm performance statistics, series limit acceleration algorithm, etc are all plugable.
For the less brave, I would appreciate any comments you might have on a quadrature library in general.
For the braver, the code can be found in the boost vault at http://www.boost-consulting.com/vault/index.php?direction=0&ord er=&directory=Math%20-%20Numerics&
Stepping in bravely (where angels fear to tread? ;-) ) At a very quick slightly-informed glance, this looks a very useful addition to the Boost numeric repertoire, espcially as it makes modern use of C++, including lambda - FORTRAN users won't like it ;-) Paul --- Paul A Bristow Prizet Farmhouse, Kendal, Cumbria UK LA8 8AB +44 1539561830 & SMS, Mobile +44 7714 330204 & SMS pbristow@hetp.u-net.com
Stepping in bravely (where angels fear to tread? ;-) )
But not rushing in, I hope...
At a very quick slightly-informed glance, this looks a very useful addition to the Boost numeric repertoire,
Sounds like encouragement to me - I shall persevere.
espcially as it makes modern use of C++, including lambda - FORTRAN users won't like it ;-)
Having been forced to read FORTRAN, I have some empathy with them :-) One of the questions I pondered was whether I should use Boost.Parameter or not. In the end I opted to use it, and I think it works well. Using thirteen parameters does however slow the compiler quite significantly. There is also a slight clash with the Math toolkit library in the vault, in that it unconditionaly defines BOOST_PARAMETER_MAX_ARITY to a smaller value. BTW, just noticed that math toolkit is in SVN - will it be in the next boost release?
Hugo Duncan wrote:
One of the questions I pondered was whether I should use Boost.Parameter or not. In the end I opted to use it, and I think it works well. Using thirteen parameters does however slow the compiler quite significantly. There is also a slight clash with the Math toolkit library in the vault, in that it unconditionaly defines BOOST_PARAMETER_MAX_ARITY to a smaller value.
That define has been removed from the Math Toolkit code in SVN.
BTW, just noticed that math toolkit is in SVN - will it be in the next boost release?
Yep :-) John.
-----Original Message----- From: boost-bounces@lists.boost.org [mailto:boost-bounces@lists.boost.org] On Behalf Of Hugo Duncan Sent: 06 December 2007 20:48 To: boost@lists.boost.org Subject: Re: [boost] [rfc] Quadrature library for numerical integrationofdefinite integrals
One of the questions I pondered was whether I should use Boost.Parameter or not. In the end I opted to use it, and I think it works well. Using thirteen parameters does however slow the compiler quite significantly.
With the GSoC Visualization using SVG written by Jake Voytko, the use of Boost.Parameter has not proved a good idea as a way of providing a zillion optional parameters. Jake is in the process of refactoring to remove it. The effect on compile time is very bad (to me entirely unacceptable) and the method of chaining (every member function returns *this) so you can write plot.color(red).title("my_title"). ... that Jake also used has proved much faster and just as convenient (in fact it allows a much nicer and more flexible layout). So I strongly recommend against using Boost.Parameter based on this experience. Paul --- Paul A Bristow Prizet Farmhouse, Kendal, Cumbria UK LA8 8AB +44 1539561830 & SMS, Mobile +44 7714 330204 & SMS pbristow@hetp.u-net.com
On Fri, 07 Dec 2007 Paul A Bristow wrote:
The effect on compile time is very bad (to me entirely unacceptable) and the method of chaining (every member function returns *this) so you can write plot.color(red).title("my_title"). ...
that Jake also used has proved much faster and just as convenient (in fact it allows a much nicer and more flexible layout).
So I strongly recommend against using Boost.Parameter based on this experience.
Heeded. Chaining could also be used to provide information on discontinuities, etc. The current functional approach does however switch between implementations at compile time depending on whether an acceleration algorithm is supplied as a parameter or not. I'll have to think about how to handle this in the chained approach ....
I have uploaded a new version that uses argument chaining. The docs have been updated to reflect this.
-----Original Message----- From: boost-bounces@lists.boost.org [mailto:boost-bounces@lists.boost.org] On Behalf Of Hugo Duncan Sent: 12 December 2007 22:22 To: boost@lists.boost.org Subject: Re: [boost] [rfc] Quadrature library for numerical integration ofdefinite integrals
I have uploaded a new version that uses argument chaining.
Even Nicer ;-) This looks a significant bit of work - definitely 'boostable' IMO ;-) It will mean we are really *starting* to get our numerical act together. Did you notice any improvement in compile time by using argument chaining? Paul PS Warning: C++ Purists may choke and splutter at finding some F**tran code in an .hpp file, but it's essential documentation from the 'bottomless pit of working subroutines'. PPS Of course it would even nicer with a templated RealType - but don't panic - version 2... --- Paul A Bristow Prizet Farmhouse, Kendal, Cumbria UK LA8 8AB +44 1539561830 & SMS, Mobile +44 7714 330204 & SMS pbristow@hetp.u-net.com PS A Typo jumped out at me std::cout << "integtral(x/log(x)) on [0,1] is " << answer
On Wed, 12 Dec 2007 Paul A Bristow wrote:
Did you notice any improvement in compile time by using argument chaining?
Not as much as I expected, though that may be due to how I have implemented the chaining - each time an argument is set, a new object is returned (with a slightly template instantiation).
PS Warning: C++ Purists may choke and splutter at finding some F**tran code in an .hpp file, but it's essential documentation from the 'bottomless pit of working subroutines'.
Oops, did I leave some in. As far as I remember, only in the integrals.hpp test file.
PPS Of course it would even nicer with a templated RealType - but don't panic - version 2...
It should already work with a templated real type. I am trying to get it to work with Boost.Interval, and have downloaded NTL to have go with that as well.
>-----Original Message----- >From: boost-bounces@lists.boost.org >[mailto:boost-bounces@lists.boost.org] On Behalf Of Hugo Duncan >Sent: 12 December 2007 23:45 >To: boost@lists.boost.org >Subject: Re: [boost] [rfc] Quadrature library for numerical >integrationofdefinite integrals > >On Wed, 12 Dec 2007 Paul A Bristow wrote: > >> Did you notice any improvement in compile time by using argument >> chaining? > >Not as much as I expected, though that may be due to how I have >implemented the chaining - each time an argument is set, a new >object is returned (with a slightly template instantiation). :-( >> PS Warning: C++ Purists may choke and splutter at finding >some F**tran >> code in an .hpp file, but it's essential documentation from >> the 'bottomless pit of working subroutines'. > >Oops, did I leave some in. As far as I remember, only in the >integrals.hpp test file. No, I was joking, but many a true word spoke in jest - it IS a vital part of the documentation. (though perhaps it need not be part of the Boost *distribution*, which reaises an interesting question about where background material should live. It needs to be somewhere *permanent* - which perhaps is not the intention or reality for the Boost Vault?) >> PPS Of course it would even nicer with a templated RealType >- but don't panic - version 2... > >It should already work with a templated real type. I am >trying to get it >to work with Boost.Interval, and have downloaded NTL to have >go with that as well. :-)) You will find essential John Maddock's patches for NTL used with the Math Toolkit - in the sandbox until 1.35 is out. And it would be really nice if we could persuade the author Victor Shoup to release at least the base part of NTL under the Boost license - is GPL at present. Good luck. Paul --- Paul A Bristow Prizet Farmhouse, Kendal, Cumbria UK LA8 8AB +44 1539561830 & SMS, Mobile +44 7714 330204 & SMS pbristow@hetp.u-net.com
Paul A Bristow wrote:
You will find essential John Maddock's patches for NTL used with the Math Toolkit - in the sandbox until 1.35 is out. And it would be really nice if we could persuade the author Victor Shoup to release at least the base part of NTL under the Boost license - is GPL at present.
Actually I've removed the patches and written a thin wrapper class instead in boost/math/bindings/rr.hpp. Hopefully folks will provide bindings for other popular packages at some point as well :-) Only trouble you might have: * NTL's error handler aborts on error: including on numeric *underflow*, something I find really annoying :-( * NTL's << streaming operator is temperamental: has a habbit of writing out the result as an integer where possible, even if the number of digits far exceeds what you can paste into C++ code as a literal. There's no way to force it into scientific mode either. Other than that it would be interesting to hear your experiences of using this. HTH, John.
Haven't studied it yet - but one question. Does this address multi-dimensional integration? Could it? Should it?
On Thu, 06 Dec 2007 Neal Becker wrote:
Haven't studied it yet - but one question. Does this address multi-dimensional integration? Could it? Should it?
It is capable of simultaneously integrating multiple functions of the same, single, variable. I use this to integrate x,y coordinates depending on an independent variable (where the cost of calculating both x and y is approximately the same as the cost of calculating either one). It does not address multi-dimensional integration. By using Fubini's theorem, it could be made to, somewhat inefficiently. Other methods (eg Monte-Carlo) have no similarities implementation wise. There are a number of other "integration" problems that could be addressed, such as ODE's, but again the algorithm's required are different. I would prefer to limit the scope to 1D integration methods. Even with 1D integrals there are areas not covered by the library - semi-infinte and infinite ranges (QUADPACK's QAGI), specification of known singularites (QAGP), Cauchy Principal Value integration (QAWC), etc
Hugo Duncan wrote:
Hi,
I have put together a quadrature library, for the numerical evaluation of definite integrals.
I put this together to meet some requirements that I had. I would like to get comments on whether the scope, algorithms and approach used are of interest within boost.
The code is based on QUADPACK's QAG and QAG global adaptive routines, using Kronrod-Gauss or Recursive Monotone Stable integration "kernels".
This is not just a translation of the QUADPACK code, but tries to take a (reasonably) modern C++ approach. The RMS kernel is based on TOMS 691.
Error estimation, choice of integration kernel, algorithm performance statistics, series limit acceleration algorithm, etc are all plugable.
For the less brave, I would appreciate any comments you might have on a quadrature library in general.
A couple of quick comments: The docs could use a quick once over with a spellchecker, I spotted "Intgeration" and "intgerators" in several places, along with "releative". It would also probably widen your audience if the docs were posted online somewhere, so folks don't have to mess about building them themselves. Other than that, can your adaptive integrators "automagically" handle functions with infinities and/or discontinuities, assuming the function is actually integratable? Or am I on a hiding to nothing expecting that? It's also good to see more Math-oriented code proposed for Boost, unfortunately I can't really claim to understand numeric integrators just yet - haven't needed them so far - but I might have a use case to test your code with if I can figure out what to do about the discontinuities! Regards, John Maddock.
The docs could use a quick once over with a spellchecker, I spotted "Intgeration" and "intgerators" in several places, along with "releative".
I'll fix these when the interface has stabilised.
It would also probably widen your audience if the docs were posted online somewhere, so folks don't have to mess about building them themselves.
The tar file contains the built docs, but I agree it would be more convenient. The docs are now available at http://www3.sympatico.ca/hugo.duncan/doc/html/
Other than that, can your adaptive integrators "automagically" handle functions with infinities and/or discontinuities, assuming the function is actually integratable? Or am I on a hiding to nothing expecting that?
It can handle some discontinuties, but not infinites (it does handle infinites at endpoints when used with kronrod_gauss). This is an area where improvement is possible.
It's also good to see more Math-oriented code proposed for Boost, unfortunately I can't really claim to understand numeric integrators just yet - haven't needed them so far - but I might have a use case to test your code with if I can figure out what to do about the discontinuities!
Are the discontinuities at known points? i.e. can you provide information to the algorithm about where the difficult points are? For my use at the moment I have functions with known discontinuities and I just split the integral at those discontinuities. Hugo
On Wed, Dec 05, 2007 at 06:00:58PM -0500, Hugo Duncan wrote:
I have put together a quadrature library, for the numerical evaluation of definite integrals.
I'm not sure whether you are interested but I would like a simple way to obtain Gaussian points for multiple dimensions for simple domains such as simplexes and of different degrees. These are the roots of the Legrende polynomials (at least in the one dimensional case) and are used for numerical integration in FEM code. Jens
On Fri, 07 Dec 2007 Jens Seidel wrote:
I'm not sure whether you are interested but I would like a simple way to obtain Gaussian points for multiple dimensions for simple domains such as simplexes and of different degrees. These are the roots of the Legrende polynomials (at least in the one dimensional case) and are used for numerical integration in FEM code.
Maybe in the future. I should have a need to do this next year. The fenics project probably has what you are loking for though http://www.fenics.org/wiki/FEniCS_Project Hugo
participants (5)
-
Hugo Duncan -
Jens Seidel -
John Maddock -
Neal Becker -
Paul A Bristow