qc25o Subroutine

public subroutine qc25o(f, a, b, omega, integr, nrmom, maxp1, ksave, result, abserr, neval, resabs, resasc, momcom, chebmo)

QC25O returns integration rules for integrands with a COS or SIN factor.

Discussion:

This routine estimates the integral I = integral of f(x) * w(x) over (a,b) where w(x) = cos(omegax) or w(x) = sin(omegax), and estimates J = integral ( A <= X <= B ) |F(X)| dx.

For small values of OMEGA or small intervals (a,b) the 15-point Gauss-Kronrod rule is used. In all other cases a generalized Clenshaw-Curtis method is used, that is, a truncated Chebyshev expansion of the function F is computed on (a,b), so that the integrand can be written as a sum of terms of the form W(X)*T(K,X), where T(K,X) is the Chebyshev polynomial of degree K. The Chebyshev moments are computed with use of a linear recurrence relation.

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, B, the limits of integration.

Input, real OMEGA, the parameter in the weight function.

Input, integer INTEGR, indicates which weight function is to be used = 1, w(x) = cos(omegax) = 2, w(x) = sin(omegax)

?, integer NRMOM, the length of interval (a,b) is equal to the length of the original integration interval divided by 2**nrmom (we suppose that the routine is used in an adaptive integration process, otherwise set nrmom = 0). nrmom must be zero at the first call.

      maxp1  - integer
               gives an upper bound on the number of Chebyshev
               moments which can be stored, i.e. for the intervals
               of lengths abs(bb-aa)*2**(-l), l = 0,1,2, ...,
               maxp1-2.

      ksave  - integer
               key which is one when the moments for the
               current interval have been computed

Output, real RESULT, the estimated value of the integral.

      abserr - real
               estimate of the modulus of the absolute
               error, which should equal or exceed abs(i-result)

Output, integer NEVAL, the number of times the integral was evaluated.

Output, real RESABS, approximation to the integral J.

Output, real RESASC, approximation to the integral of abs(F-I/(B-A)).

    on entry and return
      momcom - integer
               for each interval length we need to compute
               the Chebyshev moments. momcom counts the number
               of intervals for which these moments have already
               been computed. if nrmom < momcom or ksave = 1,
               the Chebyshev moments for the interval (a,b)
               have already been computed and stored, otherwise
               we compute them and we increase momcom.

      chebmo - real
               array of dimension at least (maxp1,25) containing
               the modified Chebyshev moments for the first momcom
               interval lengths

Local parameters:

maxp1 gives an upper bound on the number of Chebyshev moments which can be computed, i.e. for the interval (bb-aa), ..., (bb-aa)/2**(maxp1-2). should this number be altered, the first dimension of chebmo needs to be adapted.

x contains the values cos(k*pi/24) k = 1, ...,11, to be used for the Chebyshev expansion of f

      centr  - mid point of the integration interval
      hlgth  - half length of the integration interval
      fval   - value of the function f at the points
               (b-a)*0.5*cos(k*pi/12) + (b+a)*0.5
               k = 0, ...,24
      cheb12 - coefficients of the Chebyshev series expansion
               of degree 12, for the function f, in the
               interval (a,b)
      cheb24 - coefficients of the Chebyshev series expansion
               of degree 24, for the function f, in the
               interval (a,b)
      resc12 - approximation to the integral of
               cos(0.5*(b-a)*omega*x)*f(0.5*(b-a)*x+0.5*(b+a))
               over (-1,+1), using the Chebyshev series
               expansion of degree 12
      resc24 - approximation to the same integral, using the
               Chebyshev series expansion of degree 24
      ress12 - the analogue of resc12 for the sine
      ress24 - the analogue of resc24 for the sine

Arguments