Hi there. I have to evaluate an integral using GSL's quadrature algorithms, specifically, I use the QAG with 15 points. The integral has a singularity at the lower endpoint, which I have regularized following Hadamard or Gel'fand regularization technique, by subtracting 2-3 Taylor terms from the numerator of the integrand, and adding them "manually" outside the integration.
The thing is that the integrator has to be evaluated to a high precision, and thus, I couldn't use double-precision alone. Instead, I wrapped GSL with code that works with mpreal types taken from arprec, a high-precision C++ library in http://crd.lbl.gov/~dhbailey/mpdist/ . The reason for choosing arprec is because they have their own integrator, which I decided to abandon because it is simply too slow. The integrand function called by GSL's integrator essentially switches to C++ and works the integrand in high precision, where the integration variable is converted to high precision (doesn't really matter what the extra digits are, as long as the evaluation is precise), and the result is converted back to double precision before being provided to GSL. So far the settings. The questions: 1. Do you foresee and problems with the method explained? 2. For some reason, the regularization scheme doesn't work on all cases, where I get the numerator of the regularized integrand to converge to a constant non-zero value, instead of plunging down to zero faster than the denominator. Thus, the GSL integrator fails and crashes on bad behavior of the integrand, obviously. Was this observed elsewhere? Or do I sound crazy? :-) Thanx! (especially if you survived the reading thus far)! Jigal. -- Jigal A. _______________________________________________ Help-gsl mailing list [email protected] http://lists.gnu.org/mailman/listinfo/help-gsl
