[boost] Re: Re: Numerical Integration - interface

From: Dan Mcleran (dan.mcleran_at_[hidden])
Date: 2002-10-28 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
> one-dimensional 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 f-evaluation
>
> 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
trapezoidal-rule 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.

* Long-term, 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;
}
//--------------------------------------------------------------------------
-