>  Docs Center  >  Libraries  >  Markwardt  >  QPINT1D






  Craig B. Markwardt, NASA/GSFC Code 662, Greenbelt, MD 20770
  UPDATED VERSIONs can be found on my WEB PAGE:


  One dimensional numerical integration of IDL function or expression

Major Topics

  Numerical Analysis.

Calling Sequence

                  ERROR=error, NFEV=nfev, STATUS=status, NSUBINTERVALS=nsub,
                  /SYMMETRIC, SYM_AXIS= ] )


  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
  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.
  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:
      RETURN, (compute function of X and P)
  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
  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:
  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.
  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
    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
  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.


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


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


  Shows how function and expression can be used for exponential
    IDL> print, qpint1d('EXP(X)', 0D, 10D, /expr)
    IDL> print, qpint1d('EXP', 0D, 10D)
  Normal definite integral, and then parameterized using a PRIVATE
  value of 2.
    IDL> print, qpint1d('X^2*EXP(-X)', 0D, 10D, /expr)
    IDL> print, qpint1d('X^2*EXP(-X/P(0))', 0D, 10D, 2D, /expr)
  Improper integrals of the gaussian function
    IDL> inf = !values.d_infinity
    IDL> print, qpint1d('EXP(-X^2)', 0D, +inf, 2D, /expr)
    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

Common Blocks

  These common blocks are used internally only and should not be
  accessed or modified.


    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

© 2020 Harris Geospatial Solutions, Inc. |  Legal
My Account    |    Store    |    Contact Us