help-gsl
[Top][All Lists]
Advanced

[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]

[Help-gsl] Question regarding quadrature with arbitrary precision.


From: Jigal A
Subject: [Help-gsl] Question regarding quadrature with arbitrary precision.
Date: Sun, 17 Feb 2008 00:30:15 +0200

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.


reply via email to

[Prev in Thread] Current Thread [Next in Thread]