Re: [boost] Numerical Integration

From: Neal D. Becker (nbecker_at_[hidden])
Date: 2002-10-17 11:49:05

Yes I'm interested.  Here is the package I use.  You may find some
useful ideas.  I didn't propose it for boost because it uses
converted Numerical Recipes algorithms, but you could check out the
interfaces.

--------------romberg.H
#ifndef romberg_H
#define romberg_H

//! $Id: romberg.H,v 1.1 2002/04/19 12:37:06 nbecker Exp $

#include <math.h>
#include "polint.H"
#include "trapzd.H"
#include <vector>
#include <stdexcept>
#include <functional>

template<class func>
double romberg (func f, double a, double b, double eps = 1e-6, int jmax = 20, int K = 5) {
  std::vector<double> s(jmax+1);
  std::vector<double> h(jmax+2);
  
  h[1] = 1.0;
  trapzd T;
  for (int j = 1; j <= jmax; j++) {
    s[j] = T(f, a, b, j);
    if (j >= K) {
      double ss, dss;
      polint (&h[j-K], &s[j-K], K, 0.0, ss, dss);
      if (fabs (dss) <= eps * fabs (ss))
        return ss;
    }
    h[j+1] = 0.25 * h[j];
  }
  assert (1 == 0);
}

template<class func>
struct g {
  func& f;
  double ax;
  double bx;
  double eps;
  int jmax, K;

  g (func _f, double _ax, double _bx, double _eps = 1e-6, int _jmax = 20, int _K = 5) :
    f (_f), ax (_ax), bx (_bx), eps (_eps), jmax (_jmax), K(_K)
  {}

  double operator() (double y) {
    return romberg (bind2nd (f, y), ax, bx, eps, jmax, K);
  }
};

template<class func, class fa, class fb>
struct g2 {
  const func& f;
  const fa& ax;
  const fb& bx;
  double eps;
  int jmax, K;

  g2 (const func& _f, const fa& _ax, const fb& _bx, double _eps = 1e-6, int _jmax = 20, int _K = 5) :
    f (_f), ax (_ax), bx (_bx), eps (_eps), jmax (_jmax), K(_K)
  {}

  double operator() (double y) {
    return romberg (bind2nd (f, y), ax(y), bx(y), eps, jmax, K);
  }
};

template<class func>
double romberg2d (func f, double ax, double bx, double ay, double by, double eps = 1e-6, int jmax = 20, int K = 5) {
  return romberg (g<func>(f, ax, bx, eps, jmax, K), ay, by, eps, jmax, K);
}

template<class func, class fa, class fb>
double romberg2d (const func& f, const fa& ax, const fb& bx, double ay, double by, double eps = 1e-6, int jmax = 20, int K = 5) {
  return romberg (g2<func, fa, fb>(f, ax, bx, eps, jmax, K), ay, by, eps, jmax, K);
}

template<class T>
class constant : public std::unary_function<T, T> {
  T val;
public:
  constant (T _val) :
    val (_val) {}
  T operator () (const T& x) const {
    return val;
  }
};

#endif