Adaptive quadrature

This page contains the description of one-dimensional integration algorithms included into the ALGLIB package.

Adaptive quadrature: integration of smooth functions

Subroutine AutoGKSmooth includes an elementary algorithm of adaptive integration that calculates integrals of smooth functions. The algorithm is based on recursive partition of integration interval and processing the resulting subintervals with the Gauss-Kronrod quadrature formula (current implementation uses a 15-node formula). This quadrature formula allows for calculation of the integral itself and the estimate of its error. Subintervals are added into a heap and then the interval with the biggest absolute error is retrieved from the heap and partitioned again. The process is repeated until the aggregate relative error is not close to machine accuracy (the current implementation uses the following criterion: |Err| ≤ ··ε·|Integral|).

There is an alternative subroutine version - AutoGKSmoothW - which makes sure that the interval will be fragmented into subintervals that are no longer then a given value xwidth. This subroutine can be used to integrate functions with narrow peaks that can occur between two nodes of the Gauss-Kronrod quadrature formula and remain unnoticed. The subroutine AutoGKSmoothW makes sure that any peak not less then xwidth wide will be identified and taken into account when integrating the function.

Adaptive quadrature: functions with integrable singularities at the ends of the interval

The subroutines described above work well with smooth functions. However, their efficiency is substantially decreased if the function has integrable singularities at the interval ends, i.e. behaves like (x-a)^α or (b-x)^β (the singularity is integrable if α > -1 and β > -1). Sometimes thousands or tens of thousands of iterations are required to decrease the integration error to a proper level.

The AutoGKSingular subroutine solves this problem. By change of variable it can make any singularity like (x-a)^α or (b-x)^β easily integrable if the α and β coefficients are known. However, if their exact value is unknown, one can use lower bound estimate, given that the factor estimates are also strictly higher than -1.

Reverse communication and calculation of function values

The integration algorithm shall obtain values of a function during its operation. This problem is solved in most program packages by transferring the pointer to the function (C++, Delphi) or delegate (C#) which is used to calculate function value.

The ALGLIB package, differently from other libraries, makes use of reverse communication to solve this problem. When a value of a function needs to be calculated, the algorithm state is stored within a special structure, control is returned to the calling program, which makes all calculations and recalls the computing subroutine.

Thus, the integration algorithm is used in accordance in the following order:

1. AutoGKState data structure preparation by means of one of the algorithm initialization subrotuines (AutoGKSmooth, AutoGKSmoothW, AutoGKSingular).
2. AutoGKIteration subroutine call.
3. If False is returned from the subroutine, then the algorithm operation is completed and integral is successfully calculated (the integral itself can be obtained by calling the MinLMResults subprogram).
4. If True is returned from the subroutine, the latter will make a request for information on the function. Calculate function value at State.X and save it in the State.F field. Then call AutoGKIteration again.

Please note that when using the AutoGKSingular subroutine (i.e. when integrating functions with singularities at the ends of the interval) one may need the function value at the point that is too close to one of the interval ends. For instance, when integrating the f(x) = (1-x)^-0.95 function at the [0,1] interval, we may require the value f(1-10^-25). It is evident that when using double precision arithmetics, a loss of accuracy occurs, i.e. 1-10^-25 will be rounded to 1, which will lead to divizion by zero.

To avoid this situation, the following measures are taken. When it is required to calculate the function value, the algorithm fills in both State.X field and the State.XMinusA and State.BMinusX fields. In the example above, State.X and State.XMinusA will be uqual to 1 (when using a double precision arithmetics), the field State.BMinusX will have the value of 10^-25. Using the field State.BMinusX allows for avoiding the zero divide error when calculating the function value. An example of using these fields can be found in the "Subroutines Description" section.

This article is licensed for personal use only.