Re: [boost] Re: Re: Numerical Integration - interface

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