
Boost : 
From: Dan Mcleran (dan.mcleran_at_[hidden])
Date: 20021028 16:24:22
"Martin Weiser" <weiser_at_[hidden]> wrote in message
news:200210281006.53005.weiser_at_zib.de...
> On Montag, 28. Oktober 2002 00:22, Chuck Allison wrote:
> > I've been doing adaptive integration for over 20 years. It's the least
> > I think a user would expect. Maybe we should compare algorithms.
>
> Comparing algorithms is always interesting. I'd like to draw your
> attention, however, to specifying a C++ interface first. The actual
> algorithmic realization would then be a "quality of implementation"
> issue.
>
> Since there seems to be some interest in adaptive integration over
> onedimensional bounded intervals, I'll make a first stab at a possible
> interface. This is meant to be the starting point for detailed discussion
> and may contain severe errors.
>
> For the evaluation of $\int_a^b f(x) dx$, where $a,b\in\R$ and $f:\R > B$
> ($B$ should be a Banach space over the real numbers maybe even some more
> general structure) up to a certain tolerance, the following minimal
> interface could be used:
>
> template <class Functor>
> typename Functor::result_type
> integrate(typename Functor::argument_type a,
> typename Functor::argument_type b,
> Functor f,
> double tol);
>
> The following optional interface extensions may be desirable:
>  select between relative and absolute tolerance
>  specify known discontinuities
>  give a hint for the maximally useful integration order
>  select some specific algorithm
>  provide a rough estimate on flops per fevaluation
>
> Maybe this calls for a policy argument.
>
> Several theoretic properties could be made part of the interface by giving
> guarantees that the implementation will respect them at least in exact
> arithmetic (abbreviated notation):
>  linearity: integ(a,b,r*f+s*g) = r*integ(a,b,f)+s*integ(a,b,g)
>  symmetry: integ(a,b,f) = integ(b,a,f)
>  additivity: integ(a,b,f) = integ(a,c,f)+integ(c,b,f) (that's tough)
>  positivity: integ(a,b,f)>=0 if a<=b and f>=0
>
> Any thoughts?
>
> Martin
Martin,
I've expanded on your initial thoughts and written some code that performs a
trapezoidalrule integration on a couple of functors (1/x and x^2). First,
some explanation:
* I am currently not taking the error tolerance into account because I want
to see if everyone likes the interface for a 1D numerical integration class.
* Longterm, the Integrator class will be responsible for meeting the error
tolerance passed in by the user via an adaptive integration routine.
* I've expanded your interface to provide a Functor (i.e. function to
integrate) policy, and an IntegrationType policy that uses the Functor class
as a template argument.
* The IntegrationType class is responsible for integrating a single slice
under the curve by calling the Functor for a and b and then returing the
results.
* The class, Integrator, is meant to be a configurable numerical integration
class that can be configured via template arguments. It's sole job is to
ensure that the numerical integration is within the error tolerance provided
by the user.
* I wrote a couple of functors and the numerical integration class I'm
proposing inside of the hpp file.
* I used Borland C++ Builder 5, but I also compile this with g++ 2.95.2
* I have verified that the answers provided by this program match the
numerical integration utility on the web @
http://people.hofstra.edu/faculty/Stefan_Waner/RealWorld/integral/integral.h
tml
Here's the hpp file:
#ifndef _NUMERICAL_INTEGRATION_HPP_INCLUDED_
#define _NUMERICAL_INTEGRATION_HPP_INCLUDED_
#include <cmath>
//Begin Functors
template<typename T,typename U = T>class OneOver
{
public:
typedef T ParamType;
typedef U ReturnType;
public:
ReturnType operator()(const ParamType& arg)
{
return ((ReturnType)1/(ReturnType)arg);
}
};
template<typename T,typename U = T>class XSquared
{
public:
typedef T ParamType;
typedef U ReturnType;
public:
ReturnType operator()(const ParamType& arg)
{
return (arg * arg);
}
};
template<class Functor> class TrapezoidalIntegrator
{
public:
TrapezoidalIntegrator(){};
~TrapezoidalIntegrator(){};
public:
typedef typename Functor::ReturnType ReturnType;
typedef typename Functor::ParamType ParamType;
public:
ReturnType integrate(const ParamType& a,const double deltax)
{
ReturnType Ret_a = _f(a);
ReturnType Ret_b = _f((a+deltax));
double factor = (deltax/2.0);
return (factor * fabs(Ret_a + Ret_b));
}
private:
Functor _f;
};
template<class Functor,template<class>class IntegrationType> class
Integrator : public IntegrationType<Functor>
{
public:
typedef typename Functor::ParamType ParamType;
typedef typename IntegrationType<Functor>::ReturnType ReturnType;
public:
ReturnType integrate(const ParamType& a,const ParamType& b,const
double tol)
{
ReturnType area(0);
double deltax = 0.25;
for(ParamType p = a;p < b;p += deltax)
{
area += IntegrationType<Functor>::integrate(p,deltax);
}
return area;
}
};
#endif //_NUMERICAL_INTEGRATION_HPP_INCLUDED_
Here's the main.cpp (Borland Builder 5):
//

#pragma hdrstop
#include "c:\\code\\vc++\\numerical integration\\NumericalIntegration.hpp"
#include <iostream>
//

#pragma argsused
int main(int argc, char* argv[])
{
using namespace std;
unsigned char key;
double area(0.0);
Integrator<OneOver<double>,TrapezoidalIntegrator> i;
Integrator<XSquared<double>,TrapezoidalIntegrator> i2;
area = i.integrate(1.0,2.0,0.001);//integrate 1/x from 1.0 to 2.0
cout<<"Area of 1/x from 1 2 = "<<area<<endl;
area = i2.integrate(1.0,2.0,0.0001);//integrate x^2 from 1.0 to 2.0
cout<<"Area of x squared from 1  2 = "<<area<<endl;
cout<<"Press any key";
cin>>key;
return 0;
}
//

Boost list run by bdawes at acm.org, gregod at cs.rpi.edu, cpdaniel at pacbell.net, john at johnmaddock.co.uk