Maths › Calculus › Quadrature ›

Adaptive

Adaptive integration algorithm using RMS formulae

Controller: CodeCogs

Contents

Interface
Class Adaptive
Page Comments

Interface

C++

Overview

The aim of this module is to implement an algorithm that can automatically estimate the definite integral of any one-variable function over some prescribed closed interval $\inline [a, b]$ :

$I := \int_a^b f(x) \; {\rm d} x.$

(2)

The algorithm should give proper estimates for the definite integral even if the function to integrate is singular at unknown positions in the integration interval $\inline [a, b]$ .

This module implements an adaptive integration algorithm using RMS formulae, which is an improved version of the algorithm used in the quadpack Netlib library. The algorithm is able to automatically handle bad behavior in the integrand, like singularities at unknown positions in the integration interval.

Recursive monotone stable (RMS) formulae are a set of integration formulae $\inline I^{(1)}$ , $\inline I^{(2)}$ , ..., $\inline I^{(n)}$ to be applied sequentially, with the following properties:

the nodes used by the formula $\inline I^{(k - 1)}$ are a subset of the nodes used by $\inline I^{(k)}$ , for all $\inline k \in \{2, 3, \ldots, n\}$ , which means that the previously computed integral estimates are not wasted,
composition of these formulae does not imply wasting of functional values,
the precision obtained by using formula $\inline I^{(k)}$ is higher than the one obtained by using $\inline I^{(k - 1)}$ , for all $\inline k \in \{2, 3, \ldots, n\}$ ,
the distance between the nodes is shorter near the boundary of the integration interval,
all formulae are generalizations of the simple trapezoidal formula,
all the used weights are positive, i.e. all formulae are numerically stable.

The error bounds for RMS formulae are shown to be better than those for Newton-Cotes and comparable to the Clenshaw-Curtis, Gauss and Gauss-Kronrod formulae. Also, as opposed to Gauss and Gauss-Kronrod rules, RMS formulae do not miss the previously computed integral estimates.

The adaptive strategy used by this algorithm is to extrapolate the function values using the $\inline \varepsilon$ algorithm. The purpose of this algorithm is to eliminate the effects of the various kinds of singularities in the function to integrate.

This module is based on the quadpack Netlib library, but replaces its local quadrature module (LQM) from using a pair of Gauss-Kronrod formulae called QK21 , to using four RMS formulae with 13, 19, 27 and 41 nodes called $\inline Q_{13}$ , $\inline Q_{19}$ , $\inline Q_{27}$ and $\inline Q_{41}$ correspondingly. These are the best formulae found by the original authors, among the total of 27 families of RMS formulae.

The algorithm is numbered 691 in the ACM mathematical software list and its authors are Paola Favati, Grazia Lotti and Francesco Romani. A link to their paper is provided in the reference list below.

The main function in this module is integrate, which takes as parameters the integrand and the limits of the integration interval, and returns the estimated value of the integral.

The following example uses the integrate function to estimate the integral of

$f(x) = \log(x) \sqrt{\frac{x}{1 - x}}$

(3)

over the unit interval [0, 1]. Note that this integrand is singular at both end points of the unit interval.

#include <codecogs/maths/calculus/quadrature/adaptive.h>
#include <iostream>
#include <iomanip>
 
// define the integrand
double f(double x)
{
  return log(x) * sqrt(x / (1 - x));
}
 
int main()
{
  // calculate the integral
  double integral = Maths::Calculus::Quadrature::Adaptive::integrate(f, 0, 1);
 
  // set the display precision to 15 digits
  std::cout << std::setprecision(15);
 
  // display the integral value
  std::cout << "Integral: " << integral << std::endl;
 
  return 0;
}

Output:

Integral: -0.606789763508705

References

Algorithm 691: Improving QUADPACK automatic integration routines, http://portal.acm.org/citation.cfm?id=108556.108580
The quadpack library, http://www.netlib.org/quadpack/
The MathWorld page on Wynn's epsilon method, http://mathworld.wolfram.com/WynnsEpsilonMethod.html

Class Adaptive

Members of Adaptive

Integrate

static doubleintegrate(	double	(`f`)(double)[function pointer]*
	double	`a`
	double	`b`	)

This function estimates the definite integral

$I := \int_a^b f(x) \ {\rm d}x$

(4)

using an adaptive integration algorithm based on RMS formulae.

f	the function to integrate
a	the lower limit of the integration interval
b	the upper limit of the integration interval

Returns

the estimate for the definite integral of f over [a, b]

Dqxrul

static voiddqxrul(	double	(`f`)(double)[function pointer]*
	double	`xl`
	double	`xu`
	double&	`y`
	double&	`ya`
	double&	`ym`
	int&	`ke`
	int&	`k1`
	double*	`fv1`
	double*	`fv2`
	int&	`l1`
	int&	`l2`	)

This method computes

$I := \int_a^b f(x) \, {\rm d}x$

(5)

with error estimate and conditionally computes

$J := \int_a^b |f(x)| \, {\rm d}x$

(6)

by using a RMS rule.

f	function defining the integrand $\inline f(x)$
xl	lower limit of integration
xu	upper limit of integration
y	stores an approximation to the integral, computed by applying the requested RMS rule
ya	if ke = k1, then this will store an approximation to the integral $\inline J$
ym	if ke = k1, then this will store an approximation to the integral of $\inline \|f(x) - I / (\mbox{xu} - \mbox{xl})\|$ over $\inline [\mbox{xl}, \mbox{xu}]$
ke	key for choice of local integral rule; a RMS rule is used with: 13 points if ke = 2, 19 points if ke = 3, 27 points if ke = 4 and 42 points if ke = 5
k1	value of key for which the additional estimates ya and ym are to be computed
fv1	vector containing l1 saved functional values of positive abscissas
fv2	vector containing l1 saved functional values of positive abscissas
l1	number of elements in fv1
l2	number of elements in fv2

Dqxlqm

static voiddqxlqm(	double	(`f`)(double)[function pointer]*
	double	`a`
	double	`b`
	double&	`result`
	double&	`abserr`
	double&	`resabs`
	double&	`resasc`
	double*	`vr`
	double*	`vs`
	int&	`lr`
	int&	`ls`	)

This method implements the local quadrature module (LQM) used when estimating the integral. Instead of using the pair QK21 of Gauss-Kronrod rules, this LQM implements four RMS formulae of 13, 19, 27 and 41 nodes called $\inline Q_{13}$ , $\inline Q_{19}$ , $\inline Q_{27}$ and $\inline Q_{41}$ respectively.

The integrals computed by this method are

$I := \int_a^b f(x) \, {\rm d}x,$

(7)

with its error estimate and:

$J := \int_a^b |f(x)| \, {\rm d}x.$

(8)

f	function defining the integrand $\inline f(x)$
a	lower limit of integration
b	upper limit of integration
result	approximation to the integral $\inline I$
abserr	estimate of the modulus of the absolute error, which should not exceed $\inline \|I - \mbox{result}\|$
resabs	approximation to the integral $\inline J$
resasc	approximation to the integral of $\inline \|f(x) - I/(\mbox{b} - \mbox{a})\|$ over $\inline [\mbox{a}, \mbox{b}]$
vr	vector of length lr containing the saved functional values of positive abscissas
vs	vector of length ls containing the saved functional values of positive abscissas
lr	number of elements in vr
ls	number of elements in vs

Dqxrrd

static voiddqxrrd(	double	(`f`)(double)[function pointer]*
	double*	`z`
	int	`lz`
	double&	`xl`
	double&	`xu`
	double*	`r`
	double*	`s`
	int&	`lr`
	int&	`ls`	)

This method reorders the computed functional values before the bisection of an interval.

f	function defining the integrand $\inline f(x)$
z	vector containing lz saved functional values
lz	number of elements in z
xl	lower limit of interval
xu	upper limit of interval
r	vector containing lr saved functional values for the new intervals
s	vector containing ls saved functional values for the new intervals
lr	number of elements in r
ls	number of elements in s

Dqpsrt

static voiddqpsrt(	int	`last`
	int&	`maxerr`
	double&	`ermax`
	double*	`elist`
	int*	`iord`
	int	`nrmax`	)

This method maintains the descending ordering in the list of the local error estimates resulting from the interval subdivision process. At each call, two error estimates are inserted using the sequential search method, top-down for the largest error estimate and bottom-up for the smallest error estimate.

last	number of error estimates currently in the list
maxerr	points to the nrmax -th largest error estimate currently in the list
ermax	nrmax -th largest error estimate; ermax = elist [maxerr]
elist	vector of dimension last containing the error estimates
iord	vector of dimension last, whose first k elements contain pointers to the error estimates, such that elist [iord [1]], ..., elist [iord [k]] form a decreasing sequence, where k = last if last $\inline \leq$ (DIV_LIMIT / 2 + 2) and k = DIV_LIMIT + 1 - last otherwise
nrmax	value with the property that maxerr = iord [nrmax]

Dqelg

static voiddqelg(	int&	`n`
	double*	`epstab`
	double&	`result`
	double&	`abserr`
	double*	`res3la`
	int&	`nres`	)

This method determines the limit of a given sequence of approximation, by means of the $\inline \varepsilon$ algorithm of P. Wynn. An estimate of the absolute error is also given.

The condensed $\inline \varepsilon$ table is computed. Only those elements needed for the computation of the next diagonal are preserved.

n	epstab [n] contains the new element in the first column of the $\inline \varepsilon$ table
epstab	vector of dimension 52 containing the elements of the two lower diagonals of the triangular $\inline \varepsilon$ table; the elements are numbered starting at the right-hand corner of the triangle
result	resulting approximation to the integral
abserr	estimate of the absolute error computed from result and the three previous results
res3la	vector of dimension 3 containing the last three results
nres	number of calls to the routine (should be zero at first call)

Dqxgse

static voiddqxgse(	double	(`f`)(double)[function pointer]*
	double	`a`
	double	`b`
	double&	`result`
	double&	`abserr`
	int&	`ier`
	double*	`alist`
	double*	`blist`
	double*	`rlist`
	double*	`elist`
	int*	`iord`
	int&	`last`
	double*	`valp`
	double*	`valn`
	int*	`lp`
	int*	`ln`	)

This method calculates an approximation to the definite integral

$I := \int_a^b f(x) \, {\rm d}x,$

(9)

hopefully satisfying the following claim for accuracy:

$|I - \mbox{result}| \leq \max \{ \mbox{epsabs}, \mbox{epsrel} \times |I| \}.$

(10)

f	function defining the integrand $\inline f(x)$
a	lower limit of integration
b	upper limit of integration
result	approximation to the integral
abserr	estimate of the modulus of the absolute error, which should equal or exceed \| $\inline I$ - result\|
ier	if ier = 0, then the termination of this method was normal and reliable (it is assumed that the requested accuracy has been achieved); if ier > 0, there was an abnormal termination of this method, meaning that the estimates for the integral and the error are less reliable (it is assumed that the requested accuracy has not been achieved)
alist	vector of dimension at least DIV_LIMIT, the first last elements of which are the left end points of the subintervals in the partition of the given integration range [a, b]
blist	vector of dimension at least DIV_LIMIT, the first last elements of which are the right end points of the subintervals in the partition of the given integration range [a, b]
rlist	vector of dimension at least DIV_LIMIT, the first last elements of which are the integral approximations on the subintervals
elist	vector of dimension at least DIV_LIMIT, the first last elements of which are the moduli of the absolute error estimates on the subintervals
iord	vector of dimension at least DIV_LIMIT, the first k elements of which are pointers to the error estimates over the subintervals, such that elist [iord [1]], ..., elist [iord [k]] form a decreasing sequence, with k = last if last $\inline \leq$ (DIV_LIMIT / 2 + 2) and k = DIV_LIMIT + 1 - last otherwise
last	number of subintervals actually produced in the subdivision process
valp	array of dimension at least (21, DIV_LIMIT) used to save the functional values
valn	array of dimension at least (21, DIV_LIMIT) used to save the functional values
lp	vectors of dimension at least DIV_LIMIT, used to store the actual number of functional values saved in the corresponding column of valp
ln	vectors of dimension at least DIV_LIMIT, used to store the actual number of functional values saved in the corresponding column of valn