Boost logo

Boost :

From: Martin Weiser (weiser_at_[hidden])
Date: 2002-10-30 08:21:46


On Montag, 28. Oktober 2002 22:24, Dan Mcleran wrote:

> 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).

Dan,

I've written an interface framework which I think is more flexible, and
nevertheless easy to use. It's far from complete, of course, and the
implementation is clearly suboptimal, but you get the idea.

Let me know what you think.

The code below compiles and runs with GCC 3.2. No other compiler was
tested by now.

#include <iostream>
#include <iomanip>
#include <cmath>

// example integrand
struct OneOver {
  typedef double return_type;
  typedef double argument_type;

  double operator()(double x) const { return 1/x; }
};

enum tolerance_type { RELATIVE, ABSOLUTE };

// This is the first part of the interface. It stores all options for
// the integration, providing more or less sensible defaults. Options
// can be set by calling the option methods (relative_tolerance,...)
// and using or ignoring the return value. Some option methods (for
// now only "method") affect the type, hence the return value needs to
// be used.
template <class Functor, class Formula>
class integrator {
  typedef typename Functor::argument_type domain;
  typedef typename Functor::return_type range;

  

public:
  integrator(domain a, domain b, Functor f)
    : m_a(a), m_b(b), m_f(f),
    m_toltype(ABSOLUTE), m_tol(1e-3),
    m_order_hint(5) {}

  template <class NewFormula>
  integrator(integrator<Functor,NewFormula> integ):
    m_a(integ.m_a), m_b(integ.m_b), m_f(integ.m_f),
    m_toltype(integ.m_toltype), m_tol(integ.m_tol),
    m_order_hint(integ.m_order_hint) {}
  

  Formula relative_tolerance(double tol) {
    m_toltype = RELATIVE;
    m_tol = tol;
    return Formula(*this);
  }

  Formula absolute_tolerance(double tol) {
    m_toltype = ABSOLUTE;
    m_tol = tol;
    return Formula(*this);
  }

  template <template <class F> class NewFormula>
  NewFormula<Functor> method() const {
    return NewFormula<Functor>(*this);
  }

protected:
   template <class Functor2, class Formula2> friend class integrator;
  
  domain m_a;
  domain m_b;
  Functor m_f;
  tolerance_type m_toltype;
  double m_tol;
  unsigned m_order_hint;
};

  
    
  

// Sample integration method. Integration methods derive from
// integrator<> and provide an implicit conversion to the value of the
// integral. Integration methods may provide option methods on their own.
template <class Functor>
class TrapezoidalRule: public integrator<Functor,TrapezoidalRule<Functor> > {
  typedef typename Functor::argument_type domain;
  typedef typename Functor::return_type range;

public:
  template <class Formula>
  TrapezoidalRule(integrator<Functor,Formula> integ):
    integrator<Functor,TrapezoidalRule>(integ) {}
  
  operator range () const {
    range result = 0;
    
    int const n = static_cast<int>(std::ceil(std::sqrt(std::fabs(m_b-m_a)/m_tol)));
    domain xl = m_a;
    range fl = m_f(xl);
    for (int i=0; i<n; ++i) {
      domain xr = m_a + (m_b-m_a)*(i+1)/n;
      range fr = m_f(xr);

      result += (m_b-m_a)/n * (fl+fr)/2;

      xl = xr;
      fl = fr;
    }

    return result;
  }
};

// Another sample integrator, typically more accurate by a factor of
// two.
template <class Functor>
class MidpointRule: public integrator<Functor,MidpointRule<Functor> > {
  typedef typename Functor::argument_type domain;
  typedef typename Functor::return_type range;

public:
  template <class Formula>
  MidpointRule(integrator<Functor,Formula> integ):
    integrator<Functor,MidpointRule>(integ) {}
  
  operator range () const {
    range result = 0;
    
    int const n = static_cast<int>(std::ceil(std::sqrt(std::fabs(m_b-m_a)/m_tol)));
    domain xl = m_a;
    for (int i=0; i<n; ++i) {
      domain xr = m_a + (m_b-m_a)*(i+1)/n;
      range f = m_f((xl+xr)/2);

      result += (m_b-m_a)/n * f;

      xl = xr;
    }

    return result;
  }
};

// Second part of the interface, intended to be the main entry point
// for the integration library. Usage below.
template <class Functor>
TrapezoidalRule<Functor>
integrate(typename Functor::argument_type a,
          typename Functor::argument_type b,
          Functor f) {
  return TrapezoidalRule<Functor>(integrator<Functor,TrapezoidalRule<Functor> >(a,b,f));
}

// Example usage. The cast to double is necessary for disambiguating the
// call to operator <<.
int main() {
  using namespace std;
  
  cout << setprecision(12)
       << "int(1,2,1/x) = ln(2) = " << std::log(2.0) << endl
       << "ln(2)[tol 1e-2] = " << static_cast<double>(integrate(1,2,OneOver()).relative_tolerance(1e-2)) << endl
       << "ln(2)[tol 1e-3] = " << static_cast<double>(integrate(1,2,OneOver()).relative_tolerance(1e-3)) << endl
       << "ln(2)[tol 1e-4] = " << static_cast<double>(integrate(1,2,OneOver()).relative_tolerance(1e-4)) << endl
       << "ln(2)[tol 1e-5] = " << static_cast<double>(integrate(1,2,OneOver()).relative_tolerance(1e-5)) << endl
       << "ln(2)[tol 1e-6] = " << static_cast<double>(integrate(1,2,OneOver()).relative_tolerance(1e-6)) << endl
       << "ln(2)[tol 1e-7] = " << static_cast<double>(integrate(1,2,OneOver()).relative_tolerance(1e-7)) << endl
       << "ln(2)[tol 1e-8] = " << static_cast<double>(integrate(1,2,OneOver()).relative_tolerance(1e-8)) << endl;
  
  cout << setprecision(12)
       << "int(1,2,1/x) = ln(2) = " << std::log(2.0) << endl
       << "ln(2)[tol 1e-2] = " << static_cast<double>(integrate(1,2,OneOver()).relative_tolerance(1e-2).method<MidpointRule>()) << endl
       << "ln(2)[tol 1e-3] = " << static_cast<double>(integrate(1,2,OneOver()).relative_tolerance(1e-3).method<MidpointRule>()) << endl
       << "ln(2)[tol 1e-4] = " << static_cast<double>(integrate(1,2,OneOver()).relative_tolerance(1e-4).method<MidpointRule>()) << endl
       << "ln(2)[tol 1e-5] = " << static_cast<double>(integrate(1,2,OneOver()).relative_tolerance(1e-5).method<MidpointRule>()) << endl
       << "ln(2)[tol 1e-6] = " << static_cast<double>(integrate(1,2,OneOver()).relative_tolerance(1e-6).method<MidpointRule>()) << endl
       << "ln(2)[tol 1e-7] = " << static_cast<double>(integrate(1,2,OneOver()).relative_tolerance(1e-7).method<MidpointRule>()) << endl
       << "ln(2)[tol 1e-8] = " << static_cast<double>(integrate(1,2,OneOver()).relative_tolerance(1e-8).method<MidpointRule>()) << endl;
  return 0;
}

-- 
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