QDAWC

Integrates a function f(x)/(x-c) in the Cauchy principal value sense.

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)

C — Singular point. (Input)
C must not equal A or B.

RESULT — Estimate of the integral from A to B of F(X)/(X  C). (Output)

Optional Arguments

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: ERREL =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 QDAWC (F, A, B, C, RESULT [])

Specific: The specific interface names are S_QDAWC and D_QDAWC.

FORTRAN 77 Interface

Single: CALL QDAWC (F, A, B, C, ERRABS, ERRREL, RESULT, ERREST)

Double: The double precision name is DQDAWC.

Description

The routine QDAWC uses a globally adaptive scheme in an attempt to reduce the absolute error. This routine computes integrals whose integrands have the special form w(xf(x), where w(x) = 1/(x  c). If c lies in the interval of integration, then the integral is interpreted as a Cauchy principal value. A combination of modified Clenshaw-Curtis and Gauss-Kronrod formulas are employed. In addition to the general strategy described for the IMSL routine QDAG, this routine uses an extrapolation procedure known as the ɛ-algorithm. The routine QDAWC is an implementation of the subroutine QAWC, which is fully documented by Piessens et al. (1983).

Comments

1. Workspace may be explicitly provided, if desired, by use of Q2AWC/DQ2AWC. The reference is:

CALL Q2AWC (F, A, B, C, ERRABS, ERRREL, RESULT, ERREST, MAXSUB, NEVAL, NSUBIN, ALIST, BLIST, RLIST, ELIST, IORD)

The additional arguments are as follows:

MAXSUB — Number of subintervals allowed. (Input)
A value of 500 is used by QDAWC.

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 ALISTBLIST. (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.

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. If EXACT is the exact value, QDAWC attempts to find RESULT such that ABS(EXACT  RESULT.LE. MAX(ERRABSERRREL * 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 Cauchy principal value of

 

is estimated. The values of the actual and estimated error are machine dependent.

 

USE QDAWC_INT

USE UMACH_INT

 

IMPLICIT NONE

INTEGER NOUT

REAL A, ABS, ALOG, B, C, ERRABS, ERREST, ERROR, EXACT, &

F, RESULT

INTRINSIC ABS, ALOG

EXTERNAL F

! Get output unit number

CALL UMACH (2, NOUT)

! Set limits of integration and C

A = -1.0

B = 5.0

C = 0.0

! Set error tolerances

ERRABS = 0.0

CALL QDAWC (F, A, B, C, RESULT, ERRABS=ERRABS, ERREST=ERREST)

! Print results

EXACT = ALOG(125./631.)/18.

ERROR = 2*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

F = 1.0/(5.*X**3+6.0)

RETURN

END

Output

 

Computed = -0.090 Exact = -0.090

 

Error estimate = 2.022E-06 Error = 2.980E-08