Re: [boost] Numerical Quadrature Library

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?
>>>>
>>
>> How about this question?  This is perhaps the trickiest.  It's
>> tempting to
>> ignore it- but it's also solvable.
> 
> I am sorry for having ignored that point ;-) It is indeed tricky and
> I have to confess I haven't thought about this...what is your call on
> 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] = 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);
}