QPINT1D Name
QPINT1D
Author
Craig B. Markwardt, NASA/GSFC Code 662, Greenbelt, MD 20770
craigm@lheamail.gsfc.nasa.gov
UPDATED VERSIONs can be found on my WEB PAGE:
http://cow.physics.wisc.edu/~craigm/idl/idl.html
Purpose
One dimensional numerical integration of IDL function or expression
Major Topics
Numerical Analysis.
Calling Sequence
value = QPINT1D(FUNCT, A, B, [ PRIVATE, /EXPRESSION, FUNCTARGS=,
ERROR=error, NFEV=nfev, STATUS=status, NSUBINTERVALS=nsub,
EPSABS=, EPSREL=, LIMIT=, BREAKPOINTS=, NPOINTS=
/SYMMETRIC, SYM_AXIS= ] )
Description
QPINT1D adaptively calculates an approximation result to a given
definite integral
result = Integral[ f(x) dx ] over [a,b]
hopefully satisfying a constraint on the accuracy of the solution.
QPINT1D is based on the QUADPACK fortran package originally by
Piessens, de Doncker, Ueberhuber and Kahaner (and implements
equivalents to the QAGSE, QAGPE, QAGIE, and DQKxx fortran routines).
The returned result is intended to satisfy the following claim for
accuracy: ABS(result-value) LE MAX([epsabs, epsrel*ABS(value)]),
where VALUE is the true value of the integral, and EPSABS and
EPSREL are the absolute and relative error tolerances defined
below. An estimate of the error is returned in the ERROR keyword.
Either A or B may be finite or infinite (i.e., an improper
integral).
QPINT1D is "adaptive" in the sense that it locates regions of the
integration interval which contain the highest error, and focusses
its efforts on those regions. The algorithm locates these regions
by successively bisecting the starting interval. Each subinterval
is assigned an error estimate, and the region with the largest
error estimate is subdivided further, until each subinterval
carries approximately the same amount of error. Convergence of the
procedure may be accelerated by the Epsilon algorithm due to Wynn.
The estimate of the integral and the estimate of the error in each
subinterval are computed using Gauss Kronrod quadrature.
Integrators based on the 15-, 21-, 31-, 41-, 51- and 61-point
Gauss-Kronrod rule are available, and selected using the NPOINTS
keyword. Generally, the more points the greater the precision,
especially for rapidly varying functions. However the default
value of 21 is often sufficient, especially because of the adaptive
nature of QPINT1D.
In the following sections the requirements for the form of the
integrand are established. Also, a description of how QPINT1D
handles singularities and discontinuities is presented.
INTEGRAND is a FUNCTION
The integrand can be specified in two forms, either as a standard
IDL function, or as an IDL expression. If integrating a function,
then the FUNCT should be a string naming the function. The
function must be declared as following:
FUNCTION MYFUNCT, X, P, KEYWORDS=...
RETURN, (compute function of X and P)
END
The function must accept at least one, but optionally two,
parameters. The first, 'X', is a vector of abcissae where the
function is to be computed. The function must return the same
number of function values as abcissae passed. The second
positional parameter, 'P', is a purely optional PRIVATE parameter
as described below. MYFUNCT may accept more positional parameters,
but QPINT1D will not use them. The difference between X and P is
that X is the variable of integration, while P contains any other
information expected to remain essentially constant over the
integration.
INTEGRAND is an EXPRESSION
The integrand can also be specfied as an IDL expression and setting
the EXPRESSION keyword. Any expression that can accept a vector of
abcissae named 'X' and produce a corresponding vector of output is
a valid expression. Here is an example:
RESULT = QPINT1D('X^2 * EXP(-X)', /EXPRESSION, 0D, 10D)
It is important to note that the variable of integration must
always be named 'X', and the expression must be vectorizable. The
expression may also use the PRIVATE data, and as above, it would be
referred to according to the variable 'P'. For example, if the
exponential decay constant is parameterized by PRIVATE(0), then the
expression would be:
RESULT = QPINT1D('X^2 * EXP(-X/P(0))', /EXPRESSION, 0D, 10D, PRIVATE)
The user is solely responsible for defining and using the PRIVATE
data. QPINT1D does not access or modify PRIVATE / P; it only
passes it on to the user routine for convenience.
IMPROPER INTEGRALS and DISCONTINUITIES
QPINT1D computes improper integrals, as well as integrands with
discontinuities or singularities.
Improper integrals are integrals where one or both of the limits of
integration are "infinity." (Formally, these integrals are defined
by taking the limit as the integration limit tends to infinity).
QPINT1D handles a small class of such integrals, generally for
integrands that are convergent and monotonic (i.e.,
non-oscillatory, and falling off as 1/ABS(X)^2 or steeper). Such
integrals are handled by a transformation of the original interval
into the interval [0,1].
Integrals from negative infinity to positive infinity are done in
two subintervals. By default the interval is split at X EQ 0,
however this can be controlled by using the SYM_AXIS keyword.
Users should note that if the first subinterval fails the second is
not attempted, and thus the return value VALUE should not be
trusted in those cases.
Infinite integration limits are specified by using the standard
values !VALUES.F_INFINITY or !VALUES.D_INFINITY. No other special
invocation syntax is required.
The integration routine is able to handle integrands which have
integrable singularities at the endpoints. For example, the
integral:
RESULT = QPINT1D('2*sqrt((1-x)/(1+x))/(1-x^2)', 0.0d, 1d, /expr)
has a singularity at a value of X EQ 1. Still, the singularity is
integrable, and the value returned is a correct value of 2.
If known singularities are present within the interval of
integration, then users should pass the BREAKPOINTS keyword to list
the locations of these points. QPINT1D will then integrate each
subinterval separately, while still maintaining an overall error
budget.
If known discontinuities exist in the integrand, then the user may
additionally list those points using the BREAKPOINTS keyword.
It should be noted that the algorithm used is different, depending
on whether the BREAKPOINTS keyword has been specified or not (this
is the difference between the QAGSE vs. QAGPE routines in the
original FORTRAN). The algorithm *without* BREAKPOINTS is
generally thought to be more precise than *with*. Thus, it may be
worth splitting the original integration interval manually and
invoking QPINT1D without BREAKPOINTS.
Inputs
FUNCT - by default, a scalar string containing the name of an IDL
function to be integrated. See above for the formal
definition of MYFUNCT. (No default).
If the EXPRESSION keyword is set, then FUNCT is a scalar
string containing an IDL expression to be evaluated, as
described above.
A, B - a scalar number indicating the lower and upper limits of the
interval of integration (i.e., [A, B] is the interval of
integration).
PRIVATE - any optional variable to be passed on to the function to
be integrated. For functions, PRIVATE is passed as the
second positional parameter; for expressions, PRIVATE can
be referenced by the variable 'P'. QPINT1D does not
examine or alter PRIVATE. Returns
The value of the integral. If either A or B are double precision,
then the integral is computed in double precision; otherwise the
result is returned in single precision floating point.
Keyword Parameters
BREAKPOINTS - an array of numbers specifying points within the
integration interval where the integrand is
discontinuous or singular. Out of bounds points are
ignored.
Default: undefined, i.e., no such points
EPSABS - a scalar number, the absolute error tolerance requested
in the integral computation. If less than or equal to
zero, then the value is ignored.
Default: 0
EPSREL - a scalar number, the relative (i.e., fractional) error
tolerance in the integral computation. If less than or
equal to zero, then the value is ignored.
Default: 1e-4 for float, 1d-6 for double
EXPRESSION - if set, then FUNCT is an IDL expression. Otherwise,
FUNCT is an IDL function.
ERROR - upon return, this keyword contains an estimate of the
error in the computation.
FUNCTARGS - A structure which contains the parameters to be passed
to the user-supplied function specified by FUNCT via
the _EXTRA mechanism. This is the way you can pass
additional data to your user-supplied function without
using common blocks. By default, no extra parameters
are passed to the user-supplied function.
LIMIT - a scalar, the maximum number of subintervals to create
before terminating execution. Upon return, a STATUS value
of 1 indicates such an overflow occurred.
Default: 100
NFEV - upon return, this keyword contains the number of function
calls executed (i.e., the number of abcissae).
NPOINTS - a scalar, the number of Gauss Kronrod points to use in
computing the integral over a subinterval. A larger
number of points can in principle increase the precision
of the integral, but also makes the computation take
longer. Possible values are 15, 21, 31, 41, 51, and 61.
NPOINTS is rounded up to the next nearest available set,
with a maximum of 61.
Default: 21
NSUBINTERVALS - upon return, this keyword contains the number of
subintervals the integration interval was divided
into.
STATUS - upon return, the status of the integration operation is
returned in this keyword as an integer value. A value of
zero indicates success; otherwise an abnormal condition
has occurred and the returned value should be considered
erroneous or less reliable according to STATUS:
any negative number - outright failure (reserved for
future use).
-1 - the input parameters are invalid, because
epsabs LE 0 and epsrel LT max([50*EPS,0.5d-28]),
where EPS is the machine precision, or if LIMIT
is smaller than the number of BREAKPOINTS.
0 - success.
1 - maximum number of subdivisions allowed has been
achieved. One can allow more subdivisions by
increasing the value of limit (and taking the
according dimension adjustments into
account). However, if this yields no
improvement it is advised to analyze the
integrand in order to determine the integration
difficulties. If the position of a local
difficulty can be determined (i.e. singularity,
discontinuity within the interval), it should
be supplied to the routine as an element of the
vector BREAKPOINTS.
2 - The occurrence of roundoff error is detected,
which prevents the requested tolerance from
being achieved. The error may be
under-estimated.
3 - Extremely "bad" integrand behaviour occurs at
some points of the integration interval.
4 - The algorithm does not converge. Roundoff
error is detected in the extrapolation table.
It is presumed that the requested tolerance
cannot be achieved, and that the returned
result is the best which can be obtained.
5 - The integral is probably divergent, or only
slowly convergent. It must be noted that
divergence can occur with any other value of
ier GT 0.
SYM_AXIS - a scalar number, the bisection point of the real line
for improper integrals from negative infinity to
positive infinity. Otherwise ignored.
Default: 0.
Examples
Shows how function and expression can be used for exponential
integrand:
IDL> print, qpint1d('EXP(X)', 0D, 10D, /expr)
22025.466
IDL> print, qpint1d('EXP', 0D, 10D)
22025.466
Normal definite integral, and then parameterized using a PRIVATE
value of 2.
IDL> print, qpint1d('X^2*EXP(-X)', 0D, 10D, /expr)
1.9944612
IDL> print, qpint1d('X^2*EXP(-X/P(0))', 0D, 10D, 2D, /expr)
14.005568
Improper integrals of the gaussian function
IDL> inf = !values.d_infinity
IDL> print, qpint1d('EXP(-X^2)', 0D, +inf, 2D, /expr)
0.88622693
IDL> print, qpint1d('EXP(-X^2)', -inf, +inf, 2D, /expr), sqrt(!dpi)
1.7724539 1.7724539
The second integral shows the comparison to the analytic value of
SQRT(!DPI). Common Blocks
COMMON QPINT1D_MACHAR
COMMON QPINT1D_PROFILE_COMMON
COMMON QPINT1D_QKEVAL_COMMON
These common blocks are used internally only and should not be
accessed or modified. References
R. Piessens, E. deDoncker-Kapenga, C. Uberhuber, D. Kahaner
Quadpack: a Subroutine Package for Automatic Integration
Springer Verlag, 1983. Series in Computational Mathematics v.1
515.43/Q1S 100394Z
Netlib repository: http://www.netlib.org/quadpack/
SLATEC Common Mathematical Library, Version 4.1, July 1993
a comprehensive software library containing over
1400 general purpose mathematical and statistical routines
written in Fortran 77. (http://www.netlib.org/slatec/)
Modification History
Written, Feb-Jun, 2001, CM
Documented, 04 Jun, 2001, CM
Add usage message, error checking, 15 Mar 2002, CM
Correct usage message, 28 Apr 2002, CM
More error checking when user EXPRession fails, 10 Jun 2009, CM