Integrates a function of three variables with a possible internal or endpoint singularity.
F User-supplied FUNCTION to be integrated. The form is F(X, Y, Z [, ]), where
Function Return Value
F The function value. (Output)
Required Arguments
X Independent variable. (Input)
Y Independent variable. (Input)
Z Independent variable. (Input)
Optional Arguments
FCN_DATA A derived type, s_fcn_data, which may be used to pass additional information to/from the user-supplied function. For a detailed description of this argument see FCN_DATA below.
F must be declared EXTERNAL in the calling program.
A Lower limit of integration for outer dimension. (Input)
B Upper limit of integration for outer dimension. The relative values of A and B are interpreted properly. Thus if one exchanges A and B, the sign of the answer is changed. When the integrand is positive, the sign of the result is the same as the sign of B A. (Input)
G User-supplied FUNCTION to compute the lower limit of integration for the middle dimension. The form is G(X [, ]), where
Function Return Value
G The function value. (Output)
Required Arguments
X Independent variable. (Input)
Optional Arguments
FCN_DATA A derived type, s_fcn_data, which may be used to pass additional information to/from the user-supplied function. For a detailed description of this argument see FCN_DATA below.
G must be declared EXTERNAL in the calling program.
H User-supplied FUNCTION to compute the upper limit of integration for the middle dimension. The form is H(X [, ]), where
Function Return Value
H The function value. (Output)
Required Arguments
X Independent variable. (Input)
Optional Arguments
FCN_DATA A derived type, s_fcn_data, which may be used to pass additional information to/from the user-supplied function. For a detailed description of this argument see FCN_DATA below.
H must be declared EXTERNAL in the calling program
P User-supplied FUNCTION to compute the lower limit of integration for the inner dimension. The form is P(X, Y [, ]), where
Function Return Value
P The function value. (Output)
Required Arguments
X Independent variable. (Input)
Y Independent variable. (Input)
Optional Arguments
FCN_DATA A derived type, s_fcn_data, which may be used to pass additional information to/from the user-supplied function. For a detailed description of this argument see FCN_DATA below.
P must be declared EXTERNAL in the calling program.
Q User-supplied FUNCTION to compute the upper limit of integration for the inner dimension. The form is Q(X, Y [, ]), where
Function Return Value
Q The function value. (Output)
Required Arguments
X Independent variable. (Input)
Y Independent variable. (Input)
Optional Arguments
FCN_DATA A derived type, s_fcn_data, which may be used to pass additional information to/from the user-supplied function. For a detailed description of this argument see FCN_DATA below.
Q must be declared EXTERNAL in the calling program
RESULT Estimate of the integral from A to B of the integral from G(X) to H(X) of the integral from P(X,Y) to Q(X,Y) of F(X,Y,Z). (Output)
ERRABS Absolute error tolerance. See Comment 1 for a discussion of the error
tolerances. (Input)
Default: ERRABS = 0.0.
ERRFRAC A
fraction expressing the (number of correct digits of accuracy
desired)/(number of digits of achievable precision). See Comment 1 for a discussion of the error
tolerances. (Input)
Default: ERRFRAC = 0.75.
ERRREL
The error tolerance relative to the value of the integral. See Comment 1 for a
discussion of the error tolerances. (Input)
Default: ERRREL = 0.0.
ERRPOST An a
posteriori estimate of the absolute value of the error committed while
evaluating the integrand. This value may be computed during the evaluation
of the integrand. When this optional argument is used, FCN_DATA must also be
used as FCN_DATA%RDATA(1) will
be used to pass the newly calculated value of ERRPOST back from the
evaluator, F. In
this case, the user should not use FCN_DATA%RDATA(1) for
passing other data. (Input)
Default: ERRPOST = 0.0.
ERRPRIOR An a
priori estimate of the absolute value of the relative error expected to be
committed while evaluating the integrand. Changes to this value are not detected
during evaluation of the integral. (Input)
Default: ERRPRIOR = 1.19e-7 for single
precision and 2.22d-16 for double precision.
MAXFCN The
maximum number of function values to use to compute the integral.
(Input)
Default: The number of function values is not bounded.
SINGULARITY The real part of the abscissa
of a singularity or discontinuity in the innermost integrand. If this
option is used, SINGULARITY_TYPE must
also be used. (Input)
Default: It is assumed that there is no singularity in
the innermost integrand so
SINGULARITY is not
set. It is an error to set SINGULARITY without
also setting SINGULARITY_TYPE.
SINGULARITY_TYPE A signed integer
specifying the type of singularity which occurs in the innermost integrand. If
the singularity has a leading term of the form xα
where α is not an
integer, if α is large or
has the form α =
(2n-1)/2 where n is a nonnegative integer, or the singularity is
well outside the interval, set SINGULARITY_TYPE to a
positive integer. Otherwise, set SINGULARITY_TYPE to a
negative integer. (Input)
Default: It is assumed that there is no
singularity in the innermost integrand so SINGULARITY_TYPE is
not set. It is an error to set SINGULARITY_TYPE
without also setting SINGULARITY.
FCN_DATA A derived type, s_fcn_data, which may be used to pass additional information to/from the user-supplied function. The derived type, s_fcn_data, is defined as:
type s_fcn_data
real(kind(1e0)), pointer, dimension(:) :: rdata
integer, pointer, dimension(:) :: idata
end type
in module mp_types. The double precision counterpart to s_fcn_data is named d_fcn_data. The user must include a use mp_types statement in the calling program to define this derived type. (Input/Output)
NEVAL Number of function evaluations used to calculate the integral. (Output)
ERREST An estimate of the upper bound of the magnitude of the difference between RESULT and the true value of the integral. (Output)
ISTATUS A
status flag indicating the error criteria which was satisfied on exit.
ISTATUS = -1
indicates normal termination with either the absolute or relative error
tolerance criteria satisfied.
ISTATUS = -2 indicates
normal termination with neither the absolute nor the relative error tolerance
criteria satisfied, but the error tolerance based on the locally achievable
precision is satisfied.
ISTATUS = -3 indicates
normal termination with none of the error tolerance criteria satisfied.
ISTATUS = any value
other than the above indicates abnormal termination due to an error
condition. (Output)
Generic: CALL QDAG3D (F, A, B, G, H, P, Q, RESULT [, ])
Specific: The specific interface names are S_QDAG3D and D_QDAG3D.
QDAG3D, based on the JPL Library routine SINTM, approximates an iterated three-dimensional integral of the form
The integral over three dimensions is computed by repeated integration over one dimension. The integration over one dimension is estimated using quadrature formulae due to T. N. L. Patterson (1968). Patterson described a family of formulae in which the kth formula used all the integrand values used in the k-1st formula, and added 2k-1 new integrand values in an optimal way. The first formula is the midpoint rule, the second is the three point Gauss formula, and the third is the seven point Kronrod formula. Formulae of this family of higher degree had not previously been described. This program uses formulae up to k = 8.
An error estimate is obtained by comparing the values of the integral estimated by two adjacent formulae, examining differences up to the fifteenth order, integrating round-off error, integrating error declared to have been committed during computation of the integrand, integrating a first order estimate of the effect round-off error in the abscissa has on integrand values, and including errors in the limits. The latter four methods are also used to derive a bound on the achievable precision.
If the integral over an interval cannot be estimated with sufficient accuracy, the interval is subdivided. The difference table is used to discover whether the integral is difficult to compute because the integrand is too complex or has singular behavior. In the former case, the estimated error, requested error tolerance, and difference table are used to choose a step size.
In the latter case, the difference table is used in a search algorithm to find the abscissa of the singular behavior. If the singular behavior is discovered on the end of an interval, a change of independent variable is applied to reduce the strength of the singularity.
The program also uses the difference table to detect nonintegrable singularities, jump discontinuities, and computational noise.
1. The
user provides the absolute error tolerance through optional argument ERRABS. Optional
argument ERRFRAC
represents the ratio of the (number of correct digits of accuracy desired)
to (number of digits of achievable precision). Optional argument ERRREL represents the
error tolerance relative to the value of the integral. The internal value for
ERRFRAC is
bounded between .5 and 1. By default, ERRABS and ERRREL are set to 0.0
and ERRFRAC is
set to .75. These default values usually provide all the accuracy that can be
obtained efficiently.
The error tolerance relative to the value of the
integral is applied globally (over the entire region of integration) rather than
locally (one step at a time). This policy provides true control of error
relative to the value of the integral when the integrand is not sign definite,
as well as when the integrand is sign definite. To apply the criterion of error
tolerance relative to the value of the integral, the value of the integral over
the entire region, estimated without refinement of the region, is used to derive
an absolute error tolerance that may be applied locally. If the preliminary
estimate of the value of the integral is significantly in error, and the least
restrictive error tolerance is relative to the value of the integral, the cost
of computing the integral will be larger than the cost of computing the integral
to the same degree of accuracy using appropriate values of either of the other
tolerance criteria. The preliminary estimate of the integral may be
significantly in error if the integrand is not sign definite or has large
variation.
2.
Optional arguments
SINGULARITY and
SINGULARITY_TYPE provide the user with a means to give the routine
information about the location and type of any known singularity of the
innermost integrand. When an integrand appears to have singular behavior at the
end of the interval, a transformation of the variable of integration is applied
to reduce the strength of the singularity. When an integrand appears to have
singular behavior inside the interval, the abscissa of the singularity is
determined as precisely as necessary, depending on the error tolerance, and the
interval is subdivided. The discovery of singular behavior and determination of
the abscissa of singular behavior are expensive. If the user knows of the
existence of a singularity, the efficiency of computation of the integral may be
improved by requesting an immediate transformation of the independent variable
or subdivision of the interval. It is recommended that the user select these
optional arguments for all singularities, even those outside [A, B]. If the
singularity has a leading term of the form xα where α is not an
integer, if α is large or
has the form α = (2n-1)/2
where n is a nonnegative integer, or the singularity is well outside the
interval, set
SINGULARITY_TYPE to a positive value. Otherwise, set SINGULARITY_TYPE to a
negative value. The meaning of large depends on the rest of the integrand and
the length of the interval. For the typical case, a value of about 2 is
considered large. For a singularity of the form xα log
x use the above rule, even if α is an integer.
For other types of singularities make a reasonable guess based on the above. If
several similar integrals are to be computed, some experimentation may be
useful.
When SINGULARITY_TYPE is
positive, a transformation of the form
T = TA + (X
TA)2 / (TB TA) is applied, where
TA is the abscissa of the singularity and TB is the end of the
interval. If TA is outside the interval, TB will be the end of the
interval farthest from TA. If TA is inside the interval, the
interval will immediately be subdivided at TA, and both parts will be
separately integrated with TB equal to each end of the original interval,
respectively. When SINGULARITY_TYPE is
negative, a transformation of the form T = TA + (X
TA)4 / (TB TA)3 is applied,
with TA and TB as above.
If the integrand has singularities
at more than one abscissa within the region, or more than one pole near the real
axis such that the real parts are within the region of integration, then the
interval should be subdivided at the abscissa of the singularities or the real
parts of the poles, and the integrals should be computed as separate problems,
with the results summed.
The value of
is estimated.
USE QDAG3D_INT
USE UMACH_INT
IMPLICIT NONE
! Declare variables
INTEGER NOUT
REAL A, B, ERREST, F, G, H, P, Q, RESULT
EXTERNAL F, G, H, P, Q
! Get output unit number
CALL UMACH (2, NOUT)
! Set limits of integration
A = 0.0
B = 1.0
! Set singularity value and type
CALL QDAG3D ( F, A, B, G, H, P, Q, RESULT, &
ERREST=ERREST)
! Print the results
WRITE(NOUT,*) 'Result = ', RESULT
WRITE(NOUT,9999) ERREST
9999 FORMAT('Error Estimate = ', 1PE9.1)
END
REAL FUNCTION F (X, Y, Z)
REAL X, Y, Z
F = 1.0 + X + Y + 2.0*Z
RETURN
END
REAL FUNCTION G (X)
REAL X
G = 0.0
RETURN
END
REAL FUNCTION H (X)
REAL X
H = 1.0 - X
RETURN
END
REAL FUNCTION P (X, Y)
REAL X, Y
P = 0.0
RETURN
END
REAL FUNCTION Q (X, Y)
REAL X, Y
Q = 1.0 X - Y
RETURN
END
RESULT = 0.333333
Error Estimate = 1.9E-07
PHONE: 713.784.3131 FAX:713.781.9260 |