Re: [boost] Numerical Quadrature Library

From: Theodore Papadopoulo (Theodore.Papadopoulo_at_[hidden])
Date: 2006-12-05 11:42:49

Franck Stauffer wrote:
> On Dec 5, 2006, at 4:43 PM, Neal Becker wrote:
>   
>> Good to bring this up.  Note that we've discussed this before (and  
>> I posted
>> sample code).
>>
>> I haven't seen your code yet, but you mention functors must derive  
>> from some
>> base class.  That's not desirable - and I doubt it's required.
>>     
Just for comparison, here is the classes I use. It's not the same 
integration rules and incomplete (I made no effort to split the interval
of integration in the library, I do it naively in the code: there is a 
hierarchical splitting algorithm which is available but
I did not had time to recode it in C++). The basic routines are 
converted from gsl which is itself translated from
quadpack. I give it here for comparison. Two points:

1) I also have a "Function class" (included but not needed), but in fact 
everything is templated on a Functor class and there is no need
to derive from this base class. It is more a concept that states what 
methods and types a functor class must have. Actually, this is quite
like your UnaryFunction

2) One specific need I had was to be able to integrate vectorial or 
matricial functions, which is why the function has a type associated to
it which is the type of what it returns. This return type also should be 
contrained to offer some services (namely assignement, +, -, 
multiplication by a scalar,
abs,... but not done in my code).

To really define a quadrature library, we should start defining the 
concept of function (even if in the end there will not be a base class 
for that).
I think we are quite in an agreement here, but others might have a 
different view.

Then, elementary quadrature rules and their methods should be defined 
(this is basically what I provide
with GaussKronod.H). Those for simpson or trapezoidal quadrature should 
be more simple (but I have not done
them, so....).

Then, there should be higher level routines specifying how the 
elementary integrators are used to integrate a function
over an extended interval (this can be regular splitting, or 
hierarchical splitting based on error estimation). At this point, your
implementation mixes this step and the previous (mine provides only the 
previous step).

 From the begining, one might also have provision for multiple dimension 
integration....

My two cent on this problem...
Oh, I almost forgot....  there is a GSL++, I discovered it after having 
written the above classes, but it might be worth a check.

    Theo.