Slowness is the main disadvantage of methods of the classic methods: rectangle method, trapezium method and Simpson method. They are rather simple to use, but if we try to use the first two methods to find an integral with relative error 10 ^-10, the execution time will be unacceptable. The Simpson method solves the problem faster, but its speed is often not enough either.

The alternative is to use Gauss-Kronrod quadrature formula. The function is integrated over the [A, B] interval using this formula, and we get an integral value and error estimate (usually overestimated). If the estimation is too big, we can divide an interval in half and integrate each of the halves using this formula. If the total error is still too big, we bisect an interval with the biggest error. The process is repeated until the desired precision is reached.

To speed up the selection of an interval with the biggest error, the algorithm stores an intervals in the heap. This data structure provides a very quick search of the maximal element. At that point, the parameter H is passed into the algorithm. This parameter defines a size of heap, and when the algorithm starts, it allocates a memory for 4·H numbers for the heap.

Note #1
The required size of the heap depends on the transmitted function, the size of the integration domain and the desired accuracy of the answer. For example, to integrate the oscillating function sin(1000·x) on the interval [1, 0] with absolute accuracy 10 ^-6 an integrated interval must be divided into 16 subintervals, and for calculating the integral of sin(10000·x) 128 subintervals are required. It may not be necessary to divide an area for simpler functions. To begin with the size of the heap in some hundreds of elements will be reasonable.

Note #2
It is important that the algorithm analyses the absolute error, not the relative one. At that point, relative accuracy, with which an answer could be found, won't exceed computer accuracy. For example, if the integral value is equal to 10 ²⁰, it'll be impossible to calculate it with absolute error 10 ^-3 in the majority of computers.

Note #3
This algorithm outputs the value of the integral in any case even though it doesn't meet its demand for accuracy. Thus, in a previous similar case, the algorithm just overfills the heap of any size, having calculated an interval with the most available accuracy, and then returns the error code.

Subroutine description

The algorithm uses the integrable function F(x) which should be defined by the programmer. The input parameters are given below:

• A, B define an integration interval [A,B].
• Eps is a desired absolute integration accuracy.
• HeapSize is the maximum number of subintervals, on which an interval could be divided (at that point, in the subroutine is created an array with size a*HeapSize numbers for integral table on these intervals).

Output parametrs (these parametrs are set even if an integral isn't calculated with adequate accuracy):

• Integral is a computed integral value.
• HeapUsed is a number of used intervals of decomposition. Could be used for the analysis of the algorithm execution.

The result:

• True - if an integral was successfully found with adequate absolute accuracy.
• False - if the specified number of intervals of decomposition wasn't enough to find an integral with adequate accuracy.

Report a bug

C#

C# 1.0 source.
autogk61.csharp.zip - Adaptive integration

C++

C++ source.
autogk61.cpp.zip - Adaptive integration
ablas.zip - optimized basic linear algebra subroutines with SSE2 support (for C++ sources only)

Delphi

Delphi source.
Can be compiled under FPC (in Delphi compatibility mode).
autogk61.delphi.zip - Adaptive integration

Visual Basic 6

Visual Basic 6 source.
autogk61.vb6.zip - Adaptive integration

Zonnon ^beta

Zonnon source.
Zonnon is an experimental language developed at ETH Zurich.
See www.zonnon.ethz.ch for more information.
autogk61.zonnon.zip - Adaptive integration