IMSL C Math Library
int_fcn_sing_3d
Integrates a function of three variables with a possible internal or endpoint singularity.
Synopsis
#include <imsl.h>
float imsl_f_int_fcn_sing_3d (float fcn(), float a, float b, float gcn(), float hcn(), float pcn(), float qcn(), …, 0)
The type double function is imsl_d_int_fcn_sing_3d.
Required Arguments
float fcn (float x, float y, float z) (Input/Output)
User-supplied function to be integrated.
Arguments
float x (Input)
Independent variable..
float y (Input)
Independent variable.
float z (Input)
Independent variable.
Return Value
The computed function value.
float a (Input)
Lower limit of integration for outer dimension.
float b (Input)
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.
float gcn (float x) (Input/Output)
User-supplied function to compute the lower limit of integration for the middle dimension.
Arguments
float x (Input)
Independent variable.
Return Value
The computed function value at the point x.
float hcn (float x) (Input/Output)
User-supplied function to compute the upper limit of integration for the middle dimension.
Arguments
float x (Input)
Independent variable..
Return Value
The computed function value.
float pcn (float x, float y) (Input/Output)
User-supplied function to compute the lower limit of integration for the inner dimension.
Arguments
float x (Input)
Independent variable.
float y (Input)
Independent variable.
Return Value
The computed function value.
float qcn (float x, float y) (Input/Output)
User-supplied function to compute the upper limit of integration for the inner dimension.
Arguments
float x (Input)
Independent variable.
float y (Input)
Independent variable.
Return Value
The computed function value.
Return Value
An estimate of
Synopsis with Optional Arguments
#include <imsl.h>
float imsl_f_int_fcn_sing_3d (float fcn(), float a, float b, float gcn(), float hcn(), float pcn(), float qcn(),
IMSL_FCN_W_DATA, float fcn(), float *err_post, void *data,
IMSL_GCN_W_DATA, float gcn(), void *data,
IMSL_HCN_W_DATA, float hcn(), void *data,
IMSL_PCN_W_DATA, float pcn(), void *data,
IMSL_QCN_W_DATA, float qcn(), void *data,
IMSL_ERR_ABS, float err_abs,
IMSL_ERR_FRAC, float err_frac,
IMSL_ERR_REL, float err_rel,
IMSL_ERR_PRIOR, float err_prior,
IMSL_MAX_EVALS, int maxfn,
IMSL_SINGULARITY, float singularity, int singularity_type,
IMSL_N_EVALS, int *n_evals,
IMSL_ERR_EST, float *err_est,
IMSL_ISTATUS, int *istatus,
0)
Optional Arguments
IMSL_FCN_W_DATA, float fcn (float x, float y, float z, float *err_post, void *data), float *err_post, void *data (Input)
float fcn (float x, float y, float z, float *err_post, void *data) (Input)
User supplied function to be integrated, which also accepts a pointer to an a posteriori estimate of the absolute value of the error committed while evaluating the integrand, and a pointer to data that is supplied by the user. See Passing Data to User-Supplied Functions in the introduction to this manual for more details.
Arguments
float x (Input)
Independent variable.
float y (Input)
Independent variable.
float z (Input)
Independent variable.
float *err_post (Output)
An a posteriori estimate of the absolute value of the error committed while evaluating the integrand This argument provides a means for the user to have fcn compute this value as output. Although this argument must appear in the argument list of fcn, it need not be referenced in the function. See Example 2 of int_fcn_sing_1d for an example of this.
void *data (Input)
A pointer to the data to be passed to the user-supplied function.
Return Value
The computed function value.
float *err_post (Input/Output)
An a posteriori estimate of the absolute value of the error committed while evaluating the integrand. On input, the user may supply this estimate and that value will be used as the estimate thereafter provided fcn does not calculate a new value. If an a posteriori estimate of the value of the error is not known, set err_post to 0.0 on input. On output, err_post will contain either the input value set by the user or the value calculated by fcn.
void *data (Input)
A pointer to the data to be passed to the user-supplied function.
IMSL_GCN_W_DATA, float gcn (float x, void *data), void *data (Input)
float gcn (float x, void *data) (Input)
User supplied function to compute the lower limit of integration for the middle dimension which also accepts a pointer to data that is supplied by the user. See Passing Data to User-Supplied Functions in the introduction to this manual for more details.
Arguments
float x (Input)
Independent variable.
void *data (Input)
A pointer to the data to be passed to the user-supplied function.
Return Value
The computed function value at the point x.
void *data (Input)
A pointer to the data to be passed to the user-supplied function.
IMSL_HCN_W_DATA, float hcn (float x, void *data), void *data (Input)
float hcn (float x, void *data) (Input)
User supplied function to compute the upper limit of integration for the middle dimension which also accepts a pointer to data that is supplied by the user. See Passing Data to User-Supplied Functions in the introduction to this manual for more details.
Arguments
float x (Input)
Independent variable.
void *data (Input)
A pointer to the data to be passed to the user-supplied function.
Return Value
The computed function value at the point x.
void *data (Input)
A pointer to the data to be passed to the user-supplied function.
IMSL_PCN_W_DATA, float pcn (float x, float y, void *data), void *data (Input)
float pcn (float x, float y, void *data) (Input)
User supplied function to compute the lower limit of integration for the inner dimension which also accepts a pointer to data that is supplied by the user. See Passing Data to User-Supplied Functions in the introduction to this manual for more details.
Arguments
float x (Input)
Independent variable.
float y (Input)
Independent variable.
void *data (Input)
A pointer to the data to be passed to the user-supplied function.
Return Value
The computed function value.
void *data (Input)
A pointer to the data to be passed to the user-supplied function.
IMSL_QCN_W_DATA, float qcn (float x, float y, void *data), void *data (Input)
float qcn (float x, float y, void *data) (Input)
User supplied function to compute the lower limit of integration for the inner dimension which also accepts a pointer to data that is supplied by the user. See Passing Data to User-Supplied Functions in the introduction to this manual for more details.
Arguments
float x (Input)
Independent variable.
float y (Input)
Independent variable.
void *data (Input)
A pointer to the data to be passed to the user-supplied function.
Return Value
The computed function value.
void *data (Input)
A pointer to the data to be passed to the user-supplied function.
IMSL_ERR_ABS, float  err_abs (Input)
Absolute error tolerance. See Remark 1 for a discussion on the error tolerances.
Default: err_abs = 0.0
IMSL_ERR_FRAC, float err_frac (Input)
A fraction expressing the (number of correct digits of accuracy desired)/(number of digits of achievable precision). See Remark 1 for a discussion on accuracy.
Default: err_frac = 0.75
IMSL_ERR_REL, float err_rel (Input)
The error tolerance relative to the value of the integral. See Remark 1 for a discussion on the error tolerances.
Default: err_rel = 0.0
IMSL_ERR_PRIOR, float err_prior (Input)
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.
Default: err_prior = imsl_f_machine(4)
IMSL_MAX_EVALS, float maxfn (Input)
The maximum number of function evaluations to use to compute the integral.
Default: The number of function values is not bounded.
IMSL_SINGULARITY, float singularity, int singularity_type (Input)
singularity is the real part of the abscissa of a singularity or discontinuity in the innermost integrand. By default it is assumed that there is no singularity in the integrand. singularity_type is a signed integer specifying the type of singularity which occurs in the 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. Also see Remark 2.
Default: It is assumed that there is no singularity in the innermost integrand so singularity and singularity_type are not set.
IMSL_N_EVALS, int *n_evals (Output)
Number of function evaluations used to calculate the integral.
IMSL_ERR_EST, float *err_est (Output)
An estimate of the upper bound of the magnitude of the difference between the value returned by imsl_f_int_fcn_sing_3d and the true value of the integral.
IMSL_ISTATUS, int *istatus (Output)
A status flag indicating the error criteria which was satisfied on exit.
istatus
Description
-3
Indicates normal termination with either the absolute or relative error tolerance criteria satisfied.
-4
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.
-5
Indicates normal termination with none of the error tolerance criteria satisfied.
Other
Any value other than the above indicates abnormal termination due to an error condition.
Description
The function imsl_f_int_fcn_sing_3d, 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.
Remarks
Remark 1
The user provides the absolute error tolerance through optional argument IMSL_ERR_ABS. Optional argument IMSL_ERR_FRAC represents the ratio of the (number of correct digits of accuracy desired) to (number of digits of achievable precision). Optional argument IMSL_ERR_REL represents the error tolerance relative to the value of the integral. The internal value for err_frac is bounded between .5 and 1. By default, err_abs and err_rel are set to 0.0 and err_frac 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.
Remark 2
Optional argument IMSL_SINGULARITY provides 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.
Example
The value of
is estimated.
 
#include <imsl.h>
#include <stdio.h>
 
float fcn (float x, float y, float z);
float gcn (float x);
float hcn (float x);
float pcn (float x, float y);
float qcn (float x, float y);
 
 
int main(){
float a=0.0, b=1.0, errest, result;
 
result = imsl_f_int_fcn_sing_3d(fcn, a, b, gcn, hcn, pcn, qcn,
IMSL_ERR_EST, &errest, 0);
printf("The approximation to the integral is %f\n", result);
printf("The estimated error is %6.1e\n", errest);
}
 
 
float fcn (float x, float y, float z) {
return 1.0 + x + y + 2.0*z;
}
 
float gcn (float x) {
return 0.0;
}
 
float hcn (float x) {
return 1.0 - x;
}
 
float pcn (float x, float y) {
return 0.0;
}
 
float qcn (float x, float y) {
return 1.0 - x - y;
}
Output
 
The approximation to the integral is 0.333333
The estimated error is 1.9e-007
 
Fatal Errors
IMSL_NONINTEGRABLE
The integrand apparently contains a nonintegrable singularity. The abscissa of the singularity is near #. The result has been set to NaN.
IMSL_MAX_FCN_EVAL_EXCEEDED_NAN
The maximum number of function evaluations allowed, “maxfn”, has been exceeded. “maxfn” is currently set to a value of #. The result has been set to NaN.
IMSL_CRITERIA_NOT_SATISFIED
The algorithm has terminated without satisfying any of the error tolerance criteria. The error estimate is #.
IMSL_STOP_USER_FCN
Request from user supplied function to stop algorithm. User flag = "#".