QAWFE computes Fourier integrals.
Discussion:
The routine calculates an approximation RESULT to a definite integral
I = integral of FCOS(OMEGAX) or FSIN(OMEGAX) over (A,+Infinity),
hopefully satisfying
|| I - RESULT || <= EPSABS.
Author:
Robert Piessens, Elise de Doncker-Kapenger, Christian Ueberhuber, David Kahaner
Reference:
Robert Piessens, Elise de Doncker-Kapenger, Christian Ueberhuber, David Kahaner, QUADPACK, a Subroutine Package for Automatic Integration, Springer Verlag, 1983
Parameters:
Input, external real F, the name of the function routine, of the form function f ( x ) real f real x which evaluates the integrand function.
Input, real A, the lower limit of integration.
Input, real OMEGA, the parameter in the weight function.
Input, integer INTEGR, indicates which weight function is used = 1 w(x) = cos(omegax) = 2 w(x) = sin(omegax)
Input, real EPSABS, the absolute accuracy requested.
Input, integer LIMLST, an upper bound on the number of cycles. LIMLST must be at least 1. In fact, if LIMLST < 3, the routine will end with IER= 6.
Input, integer LIMIT, an upper bound on the number of subintervals allowed in the partition of each cycle, limit >= 1.
maxp1 - integer
gives an upper bound on the number of
Chebyshev moments which can be stored, i.e.
for the intervals of lengths abs(b-a)*2**(-l),
l=0,1, ..., maxp1-2, maxp1 >= 1
Output, real RESULT, the estimated value of the integral.
Output, real ABSERR, an estimate of || I - RESULT ||.
Output, integer NEVAL, the number of times the integral was evaluated.
ier - ier = 0 normal and reliable termination of
the routine. it is assumed that the
requested accuracy has been achieved.
ier > 0 abnormal termination of the routine
the estimates for integral and error
are less reliable. it is assumed that
the requested accuracy has not been
achieved.
if omega /= 0
ier = 6 the input is invalid because
(integr /= 1 and integr /= 2) or
epsabs <= 0 or limlst < 3.
result, abserr, neval, lst are set
to zero.
= 7 bad integrand behavior occurs within
one or more of the cycles. location
and type of the difficulty involved
can be determined from the vector ierlst.
here lst is the number of cycles actually
needed (see below).
ierlst(k) = 1 the maximum number of
subdivisions (= limit)
has been achieved on the
k th cycle.
= 2 occurence of roundoff
error is detected and
prevents the tolerance
imposed on the k th cycle
from being acheived.
= 3 extremely bad integrand
behavior occurs at some
points of the k th cycle.
= 4 the integration procedure
over the k th cycle does
not converge (to within the
required accuracy) due to
roundoff in the
extrapolation procedure
invoked on this cycle. it
is assumed that the result
on this interval is the
best which can be obtained.
= 5 the integral over the k th
cycle is probably divergent
or slowly convergent. it
must be noted that
divergence can occur with
any other value of
ierlst(k).
= 8 maximum number of cycles allowed
has been achieved, i.e. of subintervals
(a+(k-1)c,a+kc) where
c = (2*int(abs(omega))+1)*pi/abs(omega),
for k = 1, 2, ..., lst.
one can allow more cycles by increasing
the value of limlst (and taking the
according dimension adjustments into
account).
examine the array iwork which contains
the error flags over the cycles, in order
to eventual look for local integration
difficulties.
if the position of a local difficulty can
be determined (e.g. singularity,
discontinuity within the interval)
one will probably gain from splitting
up the interval at this point and
calling appopriate integrators on the
subranges.
= 9 the extrapolation table constructed for
convergence acceleration of the series
formed by the integral contributions
over the cycles, does not converge to
within the required accuracy.
as in the case of ier = 8, it is advised
to examine the array iwork which contains
the error flags on the cycles.
if omega = 0 and integr = 1,
the integral is calculated by means of qagi
and ier = ierlst(1) (with meaning as described
for ierlst(k), k = 1).
rslst - real
vector of dimension at least limlst
rslst(k) contains the integral contribution
over the interval (a+(k-1)c,a+kc) where
c = (2*int(abs(omega))+1)*pi/abs(omega),
k = 1, 2, ..., lst.
note that, if omega = 0, rslst(1) contains
the value of the integral over (a,infinity).
erlst - real
vector of dimension at least limlst
erlst(k) contains the error estimate
corresponding with rslst(k).
ierlst - integer
vector of dimension at least limlst
ierlst(k) contains the error flag corresponding
with rslst(k). for the meaning of the local error
flags see description of output parameter ier.
lst - integer
number of subintervals needed for the integration
if omega = 0 then lst is set to 1.
alist, blist, rlist, elist - real
vector of dimension at least limit,
iord, nnlog - integer
vector of dimension at least limit, providing
space for the quantities needed in the
subdivision process of each cycle
chebmo - real
array of dimension at least (maxp1,25),
providing space for the Chebyshev moments
needed within the cycles
Local parameters:
c1, c2 - end points of subinterval (of length
cycle)
cycle - (2*int(abs(omega))+1)*pi/abs(omega)
psum - vector of dimension at least (limexp+2)
(see routine qextr)
psum contains the part of the epsilon table
which is still needed for further computations.
each element of psum is a partial sum of
the series which should sum to the value of
the integral.
errsum - sum of error estimates over the
subintervals, calculated cumulatively
epsa - absolute tolerance requested over current
subinterval
chebmo - array containing the modified Chebyshev
moments (see also routine qc25o)
Type | Intent | Optional | Attributes | Name | ||
---|---|---|---|---|---|---|
procedure(scalar_func) | :: | f | ||||
real | :: | a | ||||
real | :: | omega | ||||
integer | :: | integr | ||||
real | :: | epsabs | ||||
integer | :: | limlst | ||||
integer | :: | limit | ||||
integer | :: | maxp1 | ||||
real | :: | result | ||||
real | :: | abserr | ||||
integer | :: | neval | ||||
integer | :: | ier | ||||
real | :: | rslst(limlst) | ||||
real | :: | erlst(limlst) | ||||
integer | :: | ierlst(limlst) | ||||
integer | :: | lst | ||||
real | :: | alist(limit) | ||||
real | :: | blist(limit) | ||||
real | :: | rlist(limit) | ||||
real | :: | elist(limit) | ||||
integer | :: | iord(limit) | ||||
integer | :: | nnlog(limit) | ||||
real | :: | chebmo(maxp1,25) |