CodeCogs

Adaptive integration algorithm using RMS formulae

Interface

#include <codecogs/maths/calculus/quadrature/adaptive.h>

using namespace Maths::Calculus::Quadrature;

static double integrate(double (*f)(double), double a, double b)
static void dqxrul(double (*f)(double), double xl, double xu, double &y, double &ya, double &ym, int &ke, int &k1, double *fv1, double *fv2, int &l1, int &l2)
static void dqxlqm(double (*f)(double), double a, double b, double &result, double &abserr, double &resabs, double &resasc, double *vr, double *vs, int &lr, int &ls)
static void dqxrrd(double (*f)(double), double *z, int lz, double &xl, double &xu, double *r, double *s, int &lr, int &ls)
static void dqpsrt(int last, int &maxerr, double &ermax, double *elist, int *iord, int nrmax)
static void dqelg(int &n, double *epstab, double &result, double &abserr, double *res3la, int &nres)
static void dqxgse(double (*f)(double), 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)

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 $[a, b]$ :

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

(1)

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

Click for more details on RMS formulae

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 $Q_{13}$ , $Q_{19}$ , $Q_{27}$ and $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}}$

(2)

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

Members of Adaptive

STATIC CLASS METHOD

integrate

This function estimates the definite integral

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

(3)

using an adaptive integration algorithm based on RMS formulae.

Parameters

the function to integrate

the lower limit of the integration interval

the upper limit of the integration interval

Returns

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

STATIC CLASS METHOD

dqxrul

This method computes

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

(4)

with error estimate and conditionally computes

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

(5)

by using a RMS rule.

Parameters

function defining the integrand $f(x)$

lower limit of integration

upper limit of integration

stores an approximation to the integral, computed by applying the requested RMS rule

if ke = k1, then this will store an approximation to the integral $J$

if ke = k1, then this will store an approximation to the integral of $|f(x) - I / (\mbox{xu} - \mbox{xl})|$ over $[\mbox{xl}, \mbox{xu}]$

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

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

number of elements in fv1

number of elements in fv2

STATIC CLASS METHOD

dqxlqm

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 $Q_{13}$ , $Q_{19}$ , $Q_{27}$ and $Q_{41}$ respectively.

The integrals computed by this method are

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

(6)

with its error estimate and:

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

(7)

Parameters

function defining the integrand $f(x)$

lower limit of integration

upper limit of integration

result

approximation to the integral $I$

abserr

estimate of the modulus of the absolute error, which should not exceed $|I - \mbox{result}|$

resabs

approximation to the integral $J$

resasc

approximation to the integral of $|f(x) - I/(\mbox{b} - \mbox{a})|$ over $[\mbox{a}, \mbox{b}]$

vector of length lr containing the saved functional values of positive abscissas

vector of length ls containing the saved functional values of positive abscissas

number of elements in vr

number of elements in vs

STATIC CLASS METHOD

dqxrrd

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

Parameters

function defining the integrand $f(x)$

vector containing lz saved functional values

number of elements in z

lower limit of interval

upper limit of interval

vector containing lr saved functional values for the new intervals

vector containing ls saved functional values for the new intervals

number of elements in r

number of elements in s

STATIC CLASS METHOD

dqpsrt

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.

Parameters

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 $\leq$ (DIV_LIMIT / 2 + 2) and k = DIV_LIMIT + 1 - last otherwise

nrmax

value with the property that maxerr = iord [nrmax]

STATIC CLASS METHOD

dqelg

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

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

Parameters

epstab [n] contains the new element in the first column of the $\varepsilon$ table

epstab

vector of dimension 52 containing the elements of the two lower diagonals of the triangular $\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)

STATIC CLASS METHOD

dqxgse

This method calculates an approximation to the definite integral

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

(8)

hopefully satisfying the following claim for accuracy:

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

(9)

Parameters

function defining the integrand $f(x)$

lower limit of integration

upper limit of integration

result

approximation to the integral

abserr

estimate of the modulus of the absolute error, which should equal or exceed | $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 $\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

vectors of dimension at least DIV_LIMIT, used to store the actual number of functional values saved in the corresponding column of valp

vectors of dimension at least DIV_LIMIT, used to store the actual number of functional values saved in the corresponding column of valn