QDAWF
Computes a Fourier integral.
Required Arguments
F — User-supplied FUNCTION to be integrated. The form is
F(X), where
X − Independent variable. (Input)
F − The function value. (Output)
F must be declared EXTERNAL in the calling program.
A — Lower limit of integration. (Input)
IWEIGH — Type of weight function used. (Input)
IWEIGH | Weight |
1 | COS(OMEGA * X) |
2 | SIN(OMEGA * X) |
OMEGA — Parameter in the weight function. (Input)
RESULT — Estimate of the integral from A to infinity of F * WEIGHT. (Output)
Optional Arguments
ERRABS — Absolute accuracy desired. (Input)
Default: ERRABS = 1.e-3 for single precision and 1.d-8 for double precision.
ERREST — Estimate of the absolute value of the error. (Output)
Default: ERREST = 1.e-3 for single precision and 1.d-8 for double precision.
FORTRAN 90 Interface
Generic: CALL QDAWF (F, A, IWEIGH, OMEGA, RESULT [, …])
Specific: The specific interface names are S_QDAWF and D_QDAWF.
FORTRAN 77 Interface
Single: CALL QDAWF (F, A, IWEIGH, OMEGA, ERRABS, RESULT, ERREST)
Double: The double precision name is DQDAWF.
Description
The routine QDAWF uses a globally adaptive scheme in an attempt to reduce the absolute error. This routine computes integrals whose integrands have the special form w(x) f(x), where w(x) is either cos ωx or sin ωx. The integration interval is always semi-infinite of the form [A, ∞]. These Fourier integrals are approximated by repeated calls to the IMSL routine QDAWO followed by extrapolation. The routine QDAWF is an implementation of the subroutine QAWF, which is fully documented by Piessens et al. (1983).
Comments
1. Workspace may be explicitly provided, if desired, by use of Q2AWF/DQ2AWF. The reference is:
CALL Q2AWF (F, A, IWEIGH, OMEGA, ERRABS, RESULT, ERREST, MAXCYL, MAXSUB, MAXCBY, NEVAL, NCYCLE, RSLIST, ERLIST, IERLST, NSUBIN, WK, IWK)
The additional arguments are as follows:
MAXSUB — Maximum number of subintervals allowed. (Input)
A value of 365 is used by QDAWF.
MAXCYL — Maximum number of cycles allowed. (Input)
MAXCYL must be at least 3. QDAWF uses 50.
MAXCBY — Maximum number of Chebyshev moments allowed. (Input)
QDAWF uses 21.
NEVAL — Number of evaluations of F. (Output)
NCYCLE — Number of cycles used. (Output)
RSLIST — Array of length MAXCYL containing the contributions to the integral over the interval (A + (k − 1) * C, A + k * C), for k = 1, …, NCYCLE. (Output)
C = (2 * INT(ABS(OMEGA)) + 1) * PI/ABS(OMEGA).
ERLIST — Array of length MAXCYL containing the error estimates for the intervals defined in RSLIST. (Output)
IERLST — Array of length MAXCYL containing error flags for the intervals defined in RSLIST. (Output)
IERLST(K) | Meaning |
|---|
1 | The maximum number of subdivisions (MAXSUB) has been achieved on the k-th cycle. |
2 | Roundoff error prevents the desired accuracy from being achieved on the k-th cycle. |
3 | Extremely bad integrand behavior occurs at some points of the k-th cycle. |
4 | Integration procedure does not converge (to the desired accuracy) due to roundoff in the extrapolation procedure on the k-th cycle. It is assumed that the result on this interval is the best that can be obtained. |
5 | Integral over the k-th cycle is divergent or slowly convergent. |
NSUBIN — Number of subintervals generated. (Output)
WK — Work array of length 4 * MAXSUB + 25 * MAXCBY.
IWK — Work array of length 2 * MAXSUB.
2. Informational errors
Type | Code | Description |
|---|
3 | 1 | Bad integrand behavior occurred in one or more cycles. |
4 | 2 | Maximum number of cycles allowed has been reached. |
3 | 3 | Extrapolation table constructed for convergence acceleration of the series formed by the integral contributions of the cycles does not converge to the requested accuracy. |
3. If EXACT is the exact value, QDAWF attempts to find RESULT such that ABS(EXACT − RESULT) .LE. ERRABS.
Example
The value of
is estimated. The values of the actual and estimated error are machine dependent. Notice that F is coded to protect for the singularity at zero.
USE QDAWF_INT
USE UMACH_INT
USE CONST_INT
IMPLICIT NONE
INTEGER IWEIGH, NOUT
REAL A, ABS, ERRABS, ERREST, ERROR, EXACT, F, &
OMEGA, PI, RESULT
INTRINSIC ABS
EXTERNAL F
! Get output unit number
CALL UMACH (2, NOUT)
! Set lower limit of integration
A = 0.0
! Select weight W(X) = COS(PI*X/2)
IWEIGH = 1
PI = CONST('PI')
OMEGA = PI/2.0
! Set error tolerance
CALL QDAWF (F, A, IWEIGH, OMEGA, RESULT, ERREST=ERREST)
! Print results
EXACT = 1.0
ERROR = ABS(RESULT-EXACT)
WRITE (NOUT,99999) RESULT, EXACT, ERREST, ERROR
99999 FORMAT (' Computed =', F8.3, 13X, ' Exact =', F8.3, /, /, &
' Error estimate =', 1PE10.3, 6X, 'Error =', 1PE10.3)
END
!
REAL FUNCTION F (X)
REAL X
REAL SQRT
INTRINSIC SQRT
IF (X .GT. 0.0) THEN
F = 1.0/SQRT(X)
ELSE
F = 0.0
END IF
RETURN
END
Output
Computed = 1.000 Exact = 1.000
Error estimate = 6.267E-04 Error = 2.205E-06
