# Boost :

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