qawfe Subroutine

public subroutine qawfe(f, a, omega, integr, epsabs, limlst, limit, maxp1, result, abserr, neval, ier, rslst, erlst, ierlst, lst, alist, blist, rlist, elist, iord, nnlog, chebmo)

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)

Arguments

Type IntentOptional 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)

Contents