 # Boost :

From: Neal Becker (ndbecker2_at_[hidden])
Date: 2006-12-06 09:23:45

Franck Stauffer wrote:

>
> On Dec 6, 2006, at 2:34 PM, Neal Becker wrote:
>
>>>
>> [...]
>>>> 3.3 Limits on some dimensions are functions of other variables?
>>>>
>>
>> tempting to
>> ignore it- but it's also solvable.
>
> I am sorry for having ignored that point ;-) It is indeed tricky and
> this point and what do you think should be achieved?

I wrote some code quite some time ago - perhaps not very sophisticated by
modern standards - but shows some example of possibilities:

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