|
|
Boost : |
From: Martin Weiser (weiser_at_[hidden])
Date: 2002-10-29 07:32:29
On Montag, 28. Oktober 2002 22:24, Dan Mcleran wrote:
> "Martin Weiser" <weiser_at_[hidden]> wrote in message
> news:200210281006.53005.weiser_at_zib.de...
>
> > 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
> >
Hi Dan,
thanks for your actual coding. I'll comment step by step.
> 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.
Well, in what way the tolerance is interpreted exactly is definitely a
quality of implementation issue. With a finite sampling of the function
to be integrated, you are bound to rely on an error estimator. The kind
of error estimator defines the conditions under which the accuracy
requirement is satisfied. (Your sample implementation below fits into
this scheme as well, and is hence a perfectly valid implementation, but
could of course be called "low quality of implementation w.r.t. accuracy"
;)
> * Long-term, the Integrator class will be responsible for meeting the
> error tolerance passed in by the user via an adaptive integration
> routine.
Ok, I'll comment on the interface below.
> * 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.
That corresponds to my "interface extension option" of selecting a
specific algorithm. Should be sensibly defaulted.
> * 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.
Ok, this looks like separating adaptivity and integration formula into
different parts of the implementation. Say, a stepsize selection scheme
in the outer loop and an fixed integration formula given as policy class.
Now I wonder how e.g. Romberg quadrature might fit into the scheme?
I'm not sure this separation is general enough to support the most
interesting algorithms.
> * 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.
Hm, this sounds as if the adaptivity should be provided by the integrator
policy, too. This would allow Romberg quadrature.
On the other hand, this scheme does not correspond to your sample
implementation, where "adaptivity" is provided by the "outer loop".
> * 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/integ
>ral.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);
> }
> };
I'd suggest to adopt the SGI STL "Adaptable Unary Function" concept here,
unless there is a good reason not to. I.e. call the types return_type and
argument_type resp.
> 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));
> }
deltax should be ParamType. Or, as I would prefer, just the right hand end
point of the interval, i.e. 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);
> }
This is inefficient as many points are evaluated twice. There needs to be
a facility to store values which are used in two incident subintervals.
> return area;
> }
> };
This seems to be the main user visible interface to the integration
algorithms. Why would a class be preferable instead of a free function?
> #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;
> }
> //---------------------------------------------------------------------
>----- -
>
>
>
>
>
>
> _______________________________________________
> Unsubscribe & other changes:
> http://lists.boost.org/mailman/listinfo.cgi/boost
-- Dr. Martin Weiser Zuse Institute Berlin weiser_at_[hidden] Scientific Computing http://www.zib.de/weiser Numerical Analysis and Modelling
Boost list run by bdawes at acm.org, gregod at cs.rpi.edu, cpdaniel at pacbell.net, john at johnmaddock.co.uk