QDAGP
Integrates a function with singularity points given.
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)
B — Upper limit of integration. (Input)
POINTS — Array of length NPTS containing breakpoints in the range of integration. (Input)
Usually these are points where the integrand has singularities.
RESULT — Estimate of the integral from A to B of F. (Output)
Optional Arguments
NPTS — Number of break points given. (Input)
Default: NPTS = size (POINTS,1).
ERRABS — Absolute accuracy desired. (Input)
Default: ERRABS = 1.e-3 for single precision and 1.d-8 for double precision.
ERRREL — Relative accuracy desired. (Input)
Default: ERRREL = 1.e-3 for single precision and 1.d-8 for double precision.
ERREST — Estimate of the absolute value of the error. (Output)
FORTRAN 90 Interface
Generic: CALL QDAGP (F, A, B, POINTS, RESULT [, …])
Specific: The specific interface names are S_QDAGP and D_QDAGP.
FORTRAN 77 Interface
Single: CALL QDAGP (F, A, B, NPTS, POINTS, ERRABS, ERRREL, RESULT, ERREST)
Double: The double precision name is DQDAGP.
Description
The routine QDAGP uses a globally adaptive scheme in order to reduce the absolute error. It initially subdivides the interval [A, B] into NPTS + 1 user-supplied subintervals and uses a 21-point Gauss-Kronrod rule to estimate the integral over each subinterval. The error for each subinterval is estimated by comparison with the 10-point Gauss quadrature rule. This routine is designed to handle endpoint as well as interior singularities. In addition to the general strategy described in the IMSL routine QDAG, this routine employs an extrapolation procedure known as the ɛ-algorithm. The routine QDAGP is an implementation of the subroutine QAGP, which is fully documented by Piessens et al. (1983).
Comments
1. Workspace may be explicitly provided, if desired, by use of Q2AGP/DQ2AGP. The reference is:
CALL Q2AGP (F, A, B, NPTS, POINTS, ERRABS, ERRREL, RESULT, ERREST, MAXSUB, NEVAL, NSUBIN, ALIST, BLIST, RLIST, ELIST, IORD, LEVEL, WK, IWK)
The additional arguments are as follows:
MAXSUB — Number of subintervals allowed. (Input)
A value of 450 is used by QDAGP.
NEVAL — Number of evaluations of F. (Output)
NSUBIN — Number of subintervals generated. (Output)
ALIST — Array of length MAXSUB containing a list of the NSUBIN left endpoints. (Output)
BLIST — Array of length MAXSUB containing a list of the NSUBIN right endpoints. (Output)
RLIST — Array of length MAXSUB containing approximations to the NSUBIN integrals over the intervals defined by ALIST, BLIST. (Output)
ELIST — Array of length MAXSUB containing the error estimates of the NSUBIN values in RLIST. (Output)
IORD — Array of length MAXSUB. (Output)
Let K be NSUBIN if NSUBIN.LE.(MAXSUB/2 + 2), MAXSUB + 1 − NSUBIN otherwise. The first K locations contain pointers to the error estimates over the subintervals, such that ELIST(IORD(1)), …, ELIST(IORD(K)) form a decreasing sequence.
LEVEL — Array of length MAXSUB, containing the subdivision levels of the subinterval. (Output)
That is, if (AA, BB) is a subinterval of (P1, P2) where P1 as well as P2 is a user-provided break point or integration limit, then (AA, BB) has level L if ABS(BB ‑ AA) = ABS(P2 ‑ P1) * 2**( ‑L).
WK — Work array of length NPTS + 2.
IWK — Work array of length NPTS + 2.
2. Informational errors
Type | Code | Description |
|---|
4 | 1 | The maximum number of subintervals allowed has been reached. |
3 | 2 | Roundoff error, preventing the requested tolerance from being achieved, has been detected. |
3 | 3 | A degradation in precision has been detected. |
3 | 4 | Roundoff error in the extrapolation table, preventing the requested tolerance from being achieved, has been detected. |
4 | 5 | Integral is probably divergent or slowly convergent. |
3. If EXACT is the exact value, QDAGP attempts to find RESULT such that ABS(EXACT − RESULT).LE.MAX(ERRABS, ERRREL * ABS(EXACT)). To specify only a relative error, set ERRABS to zero. Similarly, to specify only an absolute error, set ERRREL to zero.
Example
The value of
is estimated. The values of the actual and estimated error are machine dependent. Note that this subroutine never evaluates the user-supplied function at the user-supplied breakpoints.
USE QDAGP_INT
USE UMACH_INT
IMPLICIT NONE
INTEGER NOUT, NPTS
REAL A, ABS, ALOG, B, ERRABS, ERREST, ERROR, ERRREL, &
EXACT, F, POINTS(2), RESULT, SQRT
INTRINSIC ABS, ALOG, SQRT
EXTERNAL F
! Get output unit number
CALL UMACH (2, NOUT)
! Set limits of integration
A = 0.0
B = 3.0
! Set error tolerances
ERRABS = 0.0
ERRREL = 0.01
! Set singularity parameters
NPTS = 2
POINTS(1) = 1.0
POINTS(2) = SQRT(2.0)
CALL QDAGP (F, A, B, POINTS, RESULT, ERRABS=ERRABS, ERRREL=ERRREL, &
ERREST=ERREST)
! Print results
EXACT = 61.0*ALOG(2.0) + 77.0/4.0*ALOG(7.0) - 27.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 ABS, ALOG
INTRINSIC ABS, ALOG
F = X**3*ALOG(ABS((X*X-1.0)*(X*X-2.0)))
RETURN
END
Output
Computed = 52.741 Exact = 52.741
Error estimate = 5.062E-01 Error = 6.104E-04
