Re: Help with QAGI from QUADPACK


In article <1146151250.394479.191260@xxxxxxxxxxxxxxxxxxxxxxxxxxxx>,
"delta" <buddy101-1@xxxxxxxxxxx> writes:
While testing QAGI from QUADPACK I ran across a problem that seems to
expose a problem in QAGI. I'm wondering if anyone here has seen
similar behavior and/or knows of a fix to QAGI.
The code below attemps to integrate the Normal PDF from negative
infinity to a bound. When I set the bound to something like 33 (or
less) I get reasonable results. When the bound is set to 34 or larger
the result approaches zero, without any error condition being flagged
by QAGI.
Any help is appreciated. BTW, I really want to stick with QUADPACK, so
suggesting an alternative package is not really an option for me. I'm
more interested in understanding what is going on in QAGI.

I downloaded QAGI from netlib, and I get the same results with both the
Sun Compiler on Solaris and Compaq Visual Fortran on Windows.

Thank you,
--Delta

program test
implicit none
integer :: limit, lenw
parameter (limit=500, lenw=4*limit)
integer :: i,j, neval, ier, last, iwork(limit)
real :: err,y, bound, work(lenw)
real, external :: FCN

bound = 34.0
call qagi(fcn, bound, -1, 0.0001, 0.0, y, err, neval, ier, &
limit, lenw, last, iwork, work)

write(*,*) "QAGI: ans = ", y
write(*,*) "QAGI: ier = ",ier
end program

! Normal PDF
real function FCN(x)
implicit none
real :: rsigsq, scale, sigma, mu, pi
real, intent(in) :: x
intrinsic atan, exp, sqrt
pi=4*atan(1.)
sigma = 1.0
mu = 0.0

rsigsq=1.e0/(sqrt(2.e0)*sigma)
scale=1.e0/sqrt(2.e0*pi)
FCN=(scale/sigma)*exp(-((x-mu)*rsigsq)**2 )
! write(*,*) x, FCN
RETURN
end function


what qagi does: in transforms the infinite interval to a finite one, then applies
adaptiv quadrature using high order gauss formulas in the usual manner
(comparison of rules of different order give error estimate).
now in your case the tranformation is from ]-inf,34] to ]0,1]
by
z=1/(35-x);
x=35-1/z
dx=1/z^2 dz
now oyur integral reads
int_{0 to 1} (scale/sigma)*exp(-(35-1/z)^2*rsigsq^2) /z^2 dz
and this is evaluaeted with guass nodes (all interior to ]0,1] )
you have now a boundary singularity but the function itself is extremely
small there, and near z=1 again the integran is extremely small:
the integrator "oversees" the peak z=1/35 , thats all.
The reason is your sloppy precision requirement. I trie it with
epsrel=1.0d-10
epsabs=1.0d-14
and got



Integral = .1000000000000000E+01
error estimate = .4117845237728320E-10
number of function evaluations = 525

hth
peter
.