I have forgotten
my Password

Or login with:

  • Facebookhttp://facebook.com/
  • Googlehttps://www.google.com/accounts/o8/id
  • Yahoohttps://me.yahoo.com
get GPL

Adaptive

Adaptive integration algorithm using RMS formulae
Controller: CodeCogs

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]: 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

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

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

Class Adaptive

Members of Adaptive

Integrate

 
static doubleintegratedouble(*f)(double)[function pointer]
doublea
doubleb )
This function estimates the definite integral using an adaptive integration algorithm based on RMS formulae.
fthe function to integrate
athe lower limit of the integration interval
bthe upper limit of the integration interval

Returns

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

Dqxrul

 
static voiddqxruldouble(*f)(double)[function pointer]
doublexl
doublexu
double& y
double& ya
double& ym
int& ke
int& k1
double*fv1
double*fv2
int& l1
int& l2 )
This method computes with error estimate and conditionally computes by using a RMS rule.
ffunction defining the integrand \inline f(x)
xllower limit of integration
xuupper limit of integration
ystores an approximation to the integral, computed by applying the requested RMS rule
yaif ke = k1, then this will store an approximation to the integral \inline J
ymif 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}]
kekey 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
k1value of key for which the additional estimates ya and ym are to be computed
fv1vector containing l1 saved functional values of positive abscissas
fv2vector containing l1 saved functional values of positive abscissas
l1number of elements in fv1
l2number of elements in fv2

Dqxlqm

 
static voiddqxlqmdouble(*f)(double)[function pointer]
doublea
doubleb
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 with its error estimate and:
ffunction defining the integrand \inline f(x)
alower limit of integration
bupper limit of integration
resultapproximation to the integral \inline I
abserrestimate of the modulus of the absolute error, which should not exceed \inline |I - \mbox{result}|
resabsapproximation to the integral \inline J
resascapproximation to the integral of \inline |f(x) - I/(\mbox{b} - \mbox{a})| over \inline [\mbox{a}, \mbox{b}]
vrvector of length lr containing the saved functional values of positive abscissas
vsvector of length ls containing the saved functional values of positive abscissas
lrnumber of elements in vr
lsnumber of elements in vs

Dqxrrd

 
static voiddqxrrddouble(*f)(double)[function pointer]
double*z
intlz
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.
ffunction defining the integrand \inline f(x)
zvector containing lz saved functional values
lznumber of elements in z
xllower limit of interval
xuupper limit of interval
rvector containing lr saved functional values for the new intervals
svector containing ls saved functional values for the new intervals
lrnumber of elements in r
lsnumber of elements in s

Dqpsrt

 
static voiddqpsrtintlast
int& maxerr
double& ermax
double*elist
int*iord
intnrmax )
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.
lastnumber of error estimates currently in the list
maxerrpoints to the nrmax -th largest error estimate currently in the list
ermaxnrmax -th largest error estimate; ermax = elist [maxerr]
elistvector of dimension last containing the error estimates
iordvector 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
nrmaxvalue with the property that maxerr = iord [nrmax]

Dqelg

 
static voiddqelgint& 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.
nepstab [n] contains the new element in the first column of the \inline \varepsilon table
epstabvector 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
resultresulting approximation to the integral
abserrestimate of the absolute error computed from result and the three previous results
res3lavector of dimension 3 containing the last three results
nresnumber of calls to the routine (should be zero at first call)

Dqxgse

 
static voiddqxgsedouble(*f)(double)[function pointer]
doublea
doubleb
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 hopefully satisfying the following claim for accuracy:
ffunction defining the integrand \inline f(x)
alower limit of integration
bupper limit of integration
resultapproximation to the integral
abserrestimate of the modulus of the absolute error, which should equal or exceed |\inline I - result|
ierif 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)
alistvector 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]
blistvector 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]
rlistvector of dimension at least DIV_LIMIT, the first last elements of which are the integral approximations on the subintervals
elistvector of dimension at least DIV_LIMIT, the first last elements of which are the moduli of the absolute error estimates on the subintervals
iordvector 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
lastnumber of subintervals actually produced in the subdivision process
valparray of dimension at least (21, DIV_LIMIT) used to save the functional values
valnarray of dimension at least (21, DIV_LIMIT) used to save the functional values
lpvectors of dimension at least DIV_LIMIT, used to store the actual number of functional values saved in the corresponding column of valp
lnvectors of dimension at least DIV_LIMIT, used to store the actual number of functional values saved in the corresponding column of valn